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FOREWORD 

In  recent  years,  the  development  of  various  air  cushion  vehicles  (ACV) 
is  drawing  worldwide  interest.  In  conjunction  with  the  Navy's  Sea  Control 
Ship  Program,  the  military  application  of  various  air  cushion  vehicles  is 
steadily  increasing.  This  study  was  initiated  to  provide  the  necessary 
mathematical  model  of  a selected  vehicle  for  *se  in  a training  device  for 
the  ACV. 

Under  the  study  contract,  a mathematical  model  of  an  Amphibious  Assault 
Landing  Craft,  (AALC),  JEFF-B,  developed  by  the  Bell  Aerospace  Company  was 
provided.  The  six-degree-of-freedom  equations  of  motion  of  the  craft  were 
derived  with  respect  to  the  pilot's  seat.  The  techniques  used  in  simulating 
the  craft's  engines  and  system  components  were  conventional.  The  technique 
used  in  simulating  the  cushion  dynamics  was  very  sophisticated.  The  .iiodel 
is  suitable  for  real  time  digital  computer  activated  training  simulation. 

This  report  contains  three  sections  and  two  appendices.  The  first 
section  is  intended  as  a general  introduction  to  the  mathematical  modeling 
problem  and  unique  vehicle  characteristics.  Second  section  describes  the 
subsystems  of  the  vehicle  and  their  integration  into  the  whole  simulation. 
The  third  section  develops  the  algebraic  model  used  for  each  component  of  the 
simulation  and  the  equations  of  motion.  The  appendices  describe  mathematical 
details  of  the  derivation  of  the  vehicle-generated  water  wave  elevations. 


WEI-HUNC  YIH 
Project  Engineer 
Naval  Training  Equipment  Center 
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SECTION  I 
INTRODUCTION 


GENERAL  DESCRIPTION 

The  mathematical  model  described  herein  is  intended  for  use  in  a pilot 
training  simulator.  It  attempts  to  model  the  vehicle  responses  to  pilot 
controls  and  external  environmental  inputs  such  as  wind  and  waves  as  faith- 
fully as  can  be  done  within  the  real  time  capability  of  digital  computation. 

The  model  is  based  on  the  description  of  the  JEFF  B craft  given  in  Bell  Aero- 
space’s Preliminary  Design  Summary  Report,  (PDSR)  updated  by  verbal  and 
written  information  obtained  as  the  vehicle  construction  progresses.  The 
capability  exists  of  updating  the  model  further  as  model  and  full  scale  test 
information  is  received  and  as  the  vehicle  hardware  is  altered.  However, 
this  model  is  not  intended  for  design  or  engineering  purposes.  In  many 
cases,  such  as  the  cushion  pressure  model,  the  calculations  are  based  on 
scanty  experimental  and  analytical  evidence  that  should  not  be  taken  for  more 
than  what  it  is;  namely,  the  best  estimate  we  can  make  of  the  dynamic  charac- 
teristics of  the  particular  system.  Extrapolation  of  such  derivations  and 
curve  fits  to  other  designs  is  not  recommended.  As  further  experimental 
and  theoretical  work  is  done  at  Naval  Ships  Research  and  Development 
Center  and  elsewhere,  more  refined  procedures  will  undoubtedly  become 
available  for  various  parts  of  the  system.  The  model  is  broken  down  into 
subsystems  in  such  a way  that  updates  are  readily  incorporated. 

Many  of  the  forces  acting  on  the  vehicle  are  curve  fits  to  experimental 
data  obtained  by  Bell  Aerospace  and  used  in  their  calculations  and  preliminary 
fan  characteristics,  and  control  forces.  However,  the  basic  dynamics  of  the 
vehicle,  in  particular  the  cushion  pressure  system  and  its  interaction  with 
the  water  or  ground  surface  has  been  developed  by  the  Draper  Laboratory 
specifically  for  the  purposes  of  this  six  degree  of  freedom  simulation. 

There  are  two  ways  a simulation  problem  of  this  type  may  be  approached. 
In  most  marine  vehicle  problems  the  user  is  interested  in  motions  and  ma- 
neuvers made  around  some  "normal*1  operating  mode.  This  leads  to  an 
assumption  of  a basic  operating  mode  and  expressions  for  perturbations  in 
series  form  and,  if  the  system  is  reasonably  linear,  it  may  be  linearized 
about  the  operating  point  and  higher  terms  in  the  series  expansion  dropped. 

If  on  the  contrary  large  excursions  from  a basic  operating  point  and  a highly 
nonlinear  system  are  encountered,  the  perturbation  approach  leads  to 
difficulties  both  in  obtaining  the  higher  order  terms  and  in  including  them  in 
a real  time  computational  model.  In  such  a case  a second  approach  may  be 
made  incorporating  the  basic  physics  of  the  vehicle  and  environmental  forces. 
For  this  pilot  training  simulation  large  variations  in  speed,  large  accelera- 
tions, large  sideslip,  yaw  and  apparent  wind  angles,  and  large  variations 
ir.  control  and  machinery  forces  must  be  anticipated.  Moreover  small 
changes  in  vehicle  state,  particularly  in  changing  ground  clearance  or  pitch 
angle,  result  in  highly  nonlinear  force  characteristics  with  large  coupling 
effects  between  degrees  of  freedom.  As  a result,  the  more  basic  approach 
of  modeling  the  forces  on  the  vehicle  as  a function  of  its  present  state  and 
recent  history  has  been  chosen  rather  than  a perturbation  approach. 
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The  greatest  difference  between  the  Air  Cushion  Vehicle  and  conventional 
marine  vehicle  modeling  is  in  the  description  of  the  air  support  system  itself. 
The  cushion  airflow  and  pressure  characteristics  gre«  tly  influence  the  motion 
of  the  craft  in  all  degrees  of  freedom  both  by  exerting  direct  forces  on  the 
vehicle  and  by  altering  the  machinery  speed  and  water  surface  wave  shape. 
Much  effort  has  therefore  been  devoted  to  a model  for  the  cushion  pressures 
that  is  valid  over  the  full  range  of  vehicle  states.  The  mathematical  model 
described  herein  lor  cushion  flow  is  fundamentally  a low  Mach  number  or 
incompressible  model.  It  assumes  that  variations  in  air  density  are  small 
in  the  system  so  that  pressures  are  determined  by  fan  and  escape  area 
characteristics.  The  volume  flow  rate  into  any  cushion  must  therefore  be 
equal  to  the  escape  rate  plus  the  rate  of  change  of  cushion  volume  with  time. 
Initial  calculations  were  made  including  the  effects  of  air  compressibility 
leading  to  extremely  high  frequency  pressure  oscillations  in  response  to  step 
heave  fc-rce  inputs.  The  frequency  of  these  pressure  variations  is  of  the 
order  of  ten  times  the  natural  frequencies  of  craft  motion.  Their  incorpo- 
ration in  a real  time  digital  model  is  not  feasible  due  to  their  high  speed  and 
their  effect  on  craft  motion  is  small.  Therefore,  the  compressible  flow 
model  was  dropped  and  further  effort  was  devoted  to  refinement  of  the  in- 
compressible flow  model.  It  should  be  stressed,  however,  that  both  these 
numerical  results  and  experimental  data  < l SES  craft  indicate  that  high 
frequency  pressure  oscillations  exist  and  can  have  major  impact  on  structural 
and  machinery  design.  They  are  filtered  out  in  the  present  mathematical 
model  only  because  of  the  real  time  computational  restraint  and  their  small 
effect  on  overall  craft  motion. 


The  cushion  pressure  system  sed  on  the  Bell  Aerospace  PDSR.  A 
centrifugal  fan  on  each  side  of  the  vt.  Jle  feeds  the  forward  and  aft  cushion 
compartnmnts  through  piping,  bags,  - ad  holes.  Total  flow  through  the  fans 
of  21  x 10^  cubic  feel  per  second  corresponds  to  a fan  pressure  of  133  psf 
and  a cushion  pressure  of  109  psf  with  a clearance  of  . 25  feet  between  water 
surface  and  skirt  hemline.  At  this  condition,  approximately  40%  of  the  fan 
outlet  air  is  exited  through  the  thrust  nozzles.  As  the  vehicle  approaches  the 
ground,  clearance  decreases,  cushion  flow  decreases,  pressure  increases, 
and  nozzle  flow  increases.  To  model  these  effects  computationally  system 
characteristics  have  been  pieced  together  from  the  Bell  FDSR  and  analytical 
assumptions. 


The  fan  characteristics  are  derived  from  the  Bell  carpet  plot  shown  in 
figure  1.  Since  rpm  is  a variable  in  the  simulation,  but  a 3 dimensional 
surface  fit  is  considered  wasteful  of  computer  time,  the  fans  are  basically 
fit  with  a parabola  at  2000  rpm  and  assumed  to  follow  the  pump  affinity 
laws.  Presssure  is  assumed  proportional  to  rpm  squared  and  flow  propor- 
tional to  rpm.  The  fan  curves  are  therefore  effectively  non  dimensionalized 
based  on  their  2000  rpm  performance. 


Pressure  loss  between  fan  and  cushion  is  assumed  to  be  proportional  to 
the  square  of  the  flow  rate  to  the  particular  cushion.  Flow  through  the  thrust 
nozzle  is  assumed  proportional  to  the  square  root  of  the  fan  pressure. 
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The  most  critical  part  of  the  cushion  pressure  sys»ein  is  the  determi- 
nation of  clearance  between  the  hemline  of  the  skirts  and  the  water  surface. 
These  clearances  determine  the  escape  area  for  cushion  air  and  therefore 
the  stiffness  of  the  vehicle  in  heave,  pitch,  and  roll.  As  pointed  out  in  the 
Bell  PDSR  and  confirmed  by  our  numerical  modeling,  assumption  of  rigid 
skirts  yields  heave  stiffnesses  far  greater  than  experiment.  The  mechanism 
appears  to  be  that  as  the  vehicle  approaches  the  surface,  the  increased 
pressure  raises  the  skirt  hemline.  Therefore,  the  clearances  do  not  close 
as  much  as  the  vehicle  moves  down.  However  as  the  clearance  departs  from 
the  equilibrium  value  of  . 25  feet,  the  vehicle  stiffness  increases  dramatically, 
as  shown  by  the  Bell  experimental  results  in  figure  2.  This  highly  nonlinear 
behavior  is  evidently  caused  by  limits  to  the  amount  of  vertical  skirt  motion. 
This  behavior  has  been  modeled  by  allowing  the  skirts  to  move  upwards  with 
increasing  cushion  pressure,  but  limiting  the  motion  to  six  inches. 

Clearance  between  the  hull  bottom  plating  and  the  water  surface  is 
determined  at  35  points  as  shown  in  figure  3.  The  points  interior  to  compart- 
ments are  used  for  volume  and  rate  of  change  of  volume  computations.  The 
points  around  the  compartment  periphery  are  used  for  determination  of  gap 
sixe  and  skirt  drag.  At  each  time  step  the  height  of  the  water  surface  due  to 
seaway  and  vehicle  generated  waves  is  computed  at  the  thirty  five  points. 

Over  land  the  ground  contour  is  used.  The  heights  of  vehicle  above  water, 
skirt  clearances,  and  skirt  motion  are  then  computed.  The  relationship 
between  compartment  pressure  and  escape  flow  rate  is  then  set  up  using  a 
discharge  coefficient  of  0. 6.  The  flow  input  to  the  cushion  from  the  fan  must 
correspond  to  the  escape  flow  plus  the  change  of  compartment  volume  with 
time.  If  it  does  not,  iteration  is  performed  until  the  correct  pressures  and 
flow  rates  are  achieved. 

The  fan  flow  rate  and  pressure  lead  to  a power  absorbed  by  the  fans. 

This  power  plus  that  absorbed  by  the  propeller  on  each  side  determines  the 
power  and  torque  loading  on  the  engine  at  its  present  rpm.  If  the  engines  on 
a side  are  producing  more  torque  than  the  fan  and  propeller  absorb,  the 
shafting  accelerates  during  the  next  time  step.  Engine  torque,  power,  and 
fuel  rate  characteristics  are  derived  from  figure  4 obtained  from  Bell 
Aerospace.  The  inertia  of  the  engines,  shafting,  fan  and  propeller  system 
on  each  side  is  approximately  3. 53  ft.  lbs.  sec^  as  seen  at  turbine  speed. 

The  interaction  between  water  surface  and  vehicle  is  rather  simple 
although  the  derivation  of  the  water  surface  shape  is  extremely  sophisticated. 
Pitch,  roll,  and  heave  forces  are  derived  from  the  cushion  pressures  and 
differences  in  pressures  between  cushion  compartments.  Drag,  sideforee 
and  yaw  moments  are  derived  from  the  slope  of  the  wave  surface  under  each 
compartment,  and  the  pressure  in  the  compartment.  Included  in  the  com- 
putation of  skirt  position  is  tre  amount  of  skirt  draging  in  the  water.  This 
is  used  with  an  empirical  drag  coefficient  to  derive  skirt  viscous  drag. 

The  computation  of  water  levels  beneath  the  craft  is  based  on  a super- 
position of  waves  existing  before  the  craft  arrived  (the  seaway  and  surf) 
and  the  vehicle  generated  wave  system.  The  seaway  is  made  up  of  deep  water 
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waves  traveling  North  offshore.  As  they  approach  the  beach  which  runs 
East-West,  the  waves  gradually  change  length,  height,  and  speed  as  a function 
of  local  water  depth.  When  they  reach  very  shallow  water  the  waves  gradually 
"break".  Their  amplitude  is  reduced  and  reaches  zero  at  the  shore  line. 

The  vehicle  generated  wave  system  is  made  up  of  waves  caused  by 
the  pressure  distribution  of  the  vehicle  predominantly  during  its  recent  past. 
The  wave  height  at  each  point  under  the  craft  can  be  computed  as  the  sum  of 
effects  of  pressure,  position,  and  orientation  over  the  last  fifty  time  steps. 
This  procedure  is  used  for  two  reasons  in  this  simulation  model.  First,  the 
important  output  for  craff  motion  is  the  total  height  of  the  wave  system  under 
the  craft.  The  air  escape  and  rate  of  change  of  cushion  volume  is  dependent 
on  these  heights  and  the  forces  on  the  vehicle  result  from  the  interaction 
between  surface  profile  and  cushion  pressure.  The  usual  formulation  using 
wave  drag  and  "added  mass"  loses  this  information  entirely.  The  second 
reason  is  that  radical  maneuvers  with  large  accelerations  and  yaw  rates  are 
anticipated,  causing  the  usual  perturbation  expansions  for  wave  forces  to  be 
inadequate. 

Other  forces  on  the  vehicle  are  derived  directly  from  Bell  PDSR  plots 
and  from  their  initial  horizontal  plane  simulation.  Thl~  includes  curve  fits 
for  aerodynamic  forces  and  moments,  propeller  forces  and  torque,  and 
rudder  forces. 
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SECTION  n 


COMPONENTS  OF  THE  SIMULATION 

INTRODUCTION 

This  section  maps  out  the  subsystems  comprising  the  AALC  simu- 
lation. A functional  breakdown  is  given  and  the  relationship  of  the  compo- 
nents to  the  whole  is  described. 

FUNCTIONAL  BREAKDOWN.  The  functional  organization  of  the  AALC 
simulation  follows  figure  5.  This  simulation  consists  of  essentially  five 
parts: 

1.  Effector  and  Control  Simulation 

2.  Ship  Motion  Simulation 

3.  Sensors  Simulation 

4.  Mission  Environment  Simulation 

5.  Outputs  to  displays  and  motion  actuators 

Effector  and  Control  Simulation.  The  controls  simulated  are  those  available 
to  the  pilot  for  maneuvering  and  propelling  the  ACV  while  on  cushion.  They 
are: 

1.  The  throttle  levers  determine  the  fuel  rate  to  each  turbine  ordered 
by  the  pilot.  There  are  six  of  these  levers,  one  for  each  turbine  gas  gene- 
rator. They  are  referred  to  in  the  Bell  literature  as  N1  control,  for  gas 
generator  rpm  control.  In  this  model  these  settings  are  referred  to  as 

N1DES  (1 6),  the  6 values  of  desired  gas  generator  rpm.  These  values 

are  input  to  the  engine  simulation  of  turbine  speed  and  power.  (Item  6 in 
figure  5. ) 

2.  The  power  turbine  output  shafts  are  linked  for  the  starboard  three 
turbines.  Therefore,  only  2 controls  are  required  for  turbine  shaft  out- 
put speed  regulation.  These  two  levers,  referred  to  in  the  Bell  literature 
as  N2  controls,  are  effectively  variable  governors.  They  represent  the 
maximum  speed  the  pilot  wishes  the  output  shafts  to  turn.  Of  course  if  the 
torque  loading  is  too  high  for  the  turbines  to  reach  that  speed  at  the  present 
N1  throttle  settings,  the  N2  control  is  ineffective.  If  the  N2  control  setting 
is  exceeded,  the  throttles  are  shut  down  until  the  desired  N2  is  reached.  In 
this  report  the  N2  control  is  referred  to  an  N2GOV  (1, 2),  the  desired  upper 
limits  to  starboard  and  port  turbine  power  output  shaft  speeds.  (Item  6 in 
figure  5.) 

3.  The  rudder  angles  are  controlled  by  foot  pedals.  There  are  rudders 
behind  each  of  the  propellers,  however,  they  turn  simultaneously.  The 
ordered  rudder  angle  is  referred  to  an  RCONT  in  the  simulation.  (Item  23 
in  figure  5. ) 
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4.  The  thrust  nozzles  are  controlled  by  two  cabin  controls.  The  pilots 
steering  wheel  rotates  from  -30°  to  +30°.  This  angle,  called  ASTEER  in 
the  simulation  controls  the  thrust  nozzle  angle  through  a 180  range.  The 
other  control,  called  SWITCH,  is  a Boolean  switch.  If  its  value  is  negative, 
the  nozzles  push  within  90°  of  ahead.  If  positive  the  nozzles  push  within 
90°  of  aft.  (Items  6 and  23  in  figure  5. ) 

5.  The  propeller  blade  pitch  control  is  by  levers  for  the  starboard  and 
port  propellerr , The  controls  are  referred  to  as  /3  control  in  the  Bell  liter- 
ature. In  this  report  the  propeller  pitches  commanded  by  the  pilot  are  called 
APROP  (1,2),  ordered  blade  pitch  on  the  starboard  and  port  propellers. 

(Item  6 in  figure  5. ) 

The  effectors  cause  forces  and  moments  on  the  vehicle  which  depend 
on  their  state,  their  control,  and  their  environment.  The  effectors  simu- 
lated are  as  follows: 

1.  Propellers.  The  propeller  thrust  and  torque  as  a function  of  rpm, 
pitch  angle,  and  inflow  velocity  are  curve  fits  based  on  the  Bell  preliminary 
horizontal  plane  simulation.  (Item  7 in  figure  5. ) 

2.  Rudders.  The  rudder  forces  and  moments  as  a function  of  rudder 
angle,  propeller  slipstream  velocity,  and  vehicle  velocity  are  curve  fit 
based  on  the  Bell  preliminary  horizontal  plane  simulation^  Item  24  in 
figure  5. ) 

3.  Thrust  Nozzles.  The  forces  and  moments  exerted  on  the  vehicle  by 
the  thrust  nozzles  are  derived  from  the  efflux  of  momentum  in  the  nozzle 
air  discharge. 

4.  Engines  and  fans.  The  engine  and  fan  characteristics  are  derived 
from  data  in  the  Bell  PDSR.  (Items  7 and  26  in  figure  5. ) 

Ship  Motion  Simulation.  The  mathematical  formulation  of  the  equations  of 
motion  (Item  8 in  figure  5)  of  the  vehicle  is  in  standard  form  for  ship  mo- 
tion and  maneuvering  computations.  However,  the  determination  of  the 
hydrodynamic  forces  acting  on  the  vehicle  due  to  its  state  and  motion  is 
vastly  different  from  displacement  surface  ship  analysis.  The  two  basic 
reasons  for  this  are: 

1.  The  ACV  is  supported  by  a cushion  of  air  whose  pressure  is  a function 
of  the  motion  as  well  as  fan  and  powering  system  dynamics.  Therefore, 
relatively  simple  hydrostatics  calculations  which  determine  equilibrium 
conditions  on  a displacement  vessel  are  replaced  by  a sophisticated  pump 
and  flow  problem.  This  problem  includes  variable  system  geometry  since 
both  the  seals  and  the  water  surface  deform  in  response  to  cushion  pres- 
sure changes. 
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2.  The  hydrodynamic  interaction  between  vehicle  and  water  surface 
is  via  a boundary  condition  on  pressure  rather  than  on  normal  velocity. 
Although  this  fact  simplifies  the  mathematics  of  computation  of  the  wave 
system  and  its  interaction  with  the  vehicle,  it  does  not  allow  the  use  of  the 
mathematics  and  empirical  data  developed  for  hydrodynamic  exciting  forces 
and  restoring  forces  on  displacement  surface  ships.  The  forces  due  to  the 
water  surface  and  cushion  air  system  are  therefore  derived  from  the  inter- 
action between  cushion  fan  and  flow  relations  and  the  water  surface  profile. 
The  water  surface  at  each  point  is  the  sum  of  the  wave  height  caused  by  the 
vehicle's  pressure  distribution,  the  vehicle  generated  waves,  and  the  wave 
height  present  when  the  vehicle  arrived,  the  seaway.  If  the  vehicle  is  over 
land,  the  wave  height  calculations  are  skipped  and  the  land  elevations  used. 
(Item  1 in  figure  5. ) 


The  seaway  is  composed  of  sinusoidal  components  in  deep  water 
running  North  toward  the  beach  which  runs  East -West  (Item  2 in  figure  5} 

As  the  waves  approach  the  beach  they  gradually  shorten,  slow,  steepen, 
and  Anally  break.  Values  of  the  wave  height  at  35  points  beneath  ihe  vehicle 
are  continuously  generated  to  input  to  the  ship  forces  and  motions.  The 
seaway  model  is  incidently  capable  of  generating  wave  heights  at  any  other 
point  relative  to  the  vehicle,  if  computer  generated  TV  simulation  of  so 
complicated  a picture  should  become  feasible.  At  present,  optical  techniques 
are  anticipated  for  the  visual  simulation. 

The  vehicle  generated  wave  system  is  calculated  by  integrating  the 
waves  generated  by  the  craft's  moving  pressure  distribution  over  the  last 
fifty  time  steps.  This  is  perhaps  the  most  technically  challenging  aspect 
of  the  simulation,  and  a large  part  of  this  report  is  devoted  to  its  derivation 
and  use.  The  waveheights  generated  by  this  mathematical  model  (Item  3 
in  figure  5)  are  input  to  the  wave  calculation. 

The  wave  calculation  (Item  4 in  figure  5. ) adds  the  -wave  heights  due 
io  seaway  and  motion  to  give  the  water  surface  profile  under  the  craft  as 
a function  of  time  (Kern  5 in  figure  5). 

The  computation  of  cushion  pressure  forces  on  the  vehicle  (Item  19 
in  figure  5. ) is  simple,  although  the  computation  of  the  pressures  is  com- 
plex. Given  the  pressures  in  each  cushion  and  assuming  uniform  pressure 
in  each  cushion,  the  horizontal  plane  forces  and  moments  are  given  by  the 
pressure  and  differences  in  wave  height  at  the  skirts.  The  vehicle  is  effec- 
tively climbing  water  hills.  The  analysis  is  identical  ashore,  except  that 
the  ground  does  not  deflect  in  response  to  the  vehicle  motion.  Heave,  pitch, 
and  roll  excitation  is  derived  from  the  cushion  pressures  and  planform 
areas  of  each  compartment. 

Air  Cushion  System  Simulation.  The  dynamics  of  the  fans,  thrust  nozzles, 
cushion  compartments  and  skirts  make  up  the  air  cushion  system  model.  The 
pressure  in  each  compartment  is  determined  by  finding  an  equilibrium  between 
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the  cushion-skirt  pressure  and  flow  characteristics  which  are  functions  of 
the  skirt  gaps  and  vehicle  motion,  and  the  fan  curves,  which  are  a function 
of  pressure,  flow,  and  fan  rpms.  The  fan  curves  and  ducting  are  derived 
from  the  Bell  PDSR  (Item  26  in  figure  5),  The  seal  dynamics  (Item  16  in 
figure  5)  combine  with  the  vehicle  attitude  and  wave  heights  to  produce  the 
skirt  gaps  (Item  15  in  figure  5).  The  fan  and  skirt  gap  pressure  flow  rela- 
tions are  combined  to  find  the  pressure  in  each  cushion  (Item  18  in  figure  5). 
These  pressures  and  flows  are  used  to  determine  nozzle  thrusts  and  intake 
air  momentum  drag  (Item  20  in  figure  5)  as  well  as  the  forces  due  directly 
to  the  pressure  distribution  (Item  19  in  figure  5). 

Sensors  Simulation.  Sensors  of  ship  motion  and  external  environment 
available  to  the  pilot  are  included  in  the  mathematical  model  and  output  at 
each  time  step  for  D/A  conversion.  These  sensors  include: 

1.  Heading  Angle.  Rotation  clockwise  from  North  is  integrated  from 
the  equations  of  motion  and  is  called  X(o)  in  the  simulation  (Item  9 in 
figure  5). 

2.  Ahead  speed  is  output  in  feet  per  second  from  integrations  of  the 
equations  of  motion  in  the  boat  coordinate  system.  Small  pitch  angles  are 
assumed  and  this  value  is  not  corrected  for  pitch.  It  is  called  XD(1)  in 
the  simulation  (Rem  1C  in  figure  5). 

3.  Sideslip  speed  is  output  in  feet  per  second  from  integration  of  the 
equations  of  motion  in  the  boat  coordinate  system.  Roll  angles  are  assumed 
to  be  small  and  this  value  is  not  corrected  for  roll  for  the  indicator  input. 

It  is  called  XD(2)  in  the  simulation  (Item  10  in  figure  5). 

4.  A height  above  water  transducer  is  assumed  attached  to  the  hull 
bottom  plating  at  the  center  of  the  ship.  This  measurement  is  available 
with  height  values  at  34  other  points  from  the  simulation  if  needed.  It  is 
called  HOW(8)  in  the  simulation  (Item  5 in  figure  5). 

5.  Roll  Angle.  The  angle  of  roll  as  seen  in  the  ship's  coordinate 
system  is  integrated  from  the  equations  of  motion.  It  is  called  X(4)  in  the 
simulation  (Item  9 in  figure  5). 

6.  Pitch  Anple.  The  angle  of  pitch  is  integrated  from  the  equations  of 
motion  and  is  called  X(5)  in  the  simulation  (Item  9 in  figure  5*. 

7.  The  angle  of  the  rudders  is  calculated  as  output  from  the  rudder 
control  (Item  23  in  figure  5)  and  is  called  RANGLE  in  the  simulation 

8.  The  thrust  nozzle  angles  are  also  calculated  as  output  from  direction- 
al control  (Item  23  in  figure  5)  and  are  called  PS  in  the  simulation. 
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9.  The  velocity  of  the  wind  relative  to  the  boat  is  calculated  from  the 
input  wind  and  the  craft  horizontal  plane  velocities  (Items  13  and  12  in 
figure  5).  It  is  called  VA  in  the  simulation.  At  the  same  time,  the  appar- 
ent wind  angle,  clockwise  from  the  bow,  is  output,  it  is  called  A BETA. 


10.  The  propeller  pitch  angles  are  output  from  the  propulsion  control 
logic  (Item  6 in  figure  5).  They  are  called  BETA(l)  and  BETA(2)  in 
the  simulation  for  starboard  and  port  propellers. 


11.  The  values  of  gas  generator  rpm  are  output  from  the  propulsion 
forces  calculation  (Item  7 in  figure  5).  Nl(l),  Nl(2),  and  Nl{3)  are  the 
values  for  the  starboard  side  turbines  and  Nl(4)  to  Nl(6)  are  the  port  side 
speeds.  All  engine  calculations  are  done  by  a subroutine  called  ENGINE. 


12.  The  values  for  starboard  and  port  power  turbine  rpm  are  also  calcu- 
lated in  the  propulsion  forces  section  (Item  7 in  figure  5.  ) They  are  called 
N2(l)  and  N2(2).  These  port  and  starboard  power  turbine  speeds  may  also 
be  used  to  indicate  fan  and  propeller  rpms  since  they  are  connected  by  gearing. 


13.  The  fuel  flow  rates  in  gallons  per  hour  output  from  the  propulsion 

forces  section  (Item  7 in  figure  5).  They  are  directly  related  to  Nl(l 6), 

the  gas  generator  speeds.  They  are  called  GFH(l)  to  GPH{6).  The  time 
integral  of  these  rates  gives  total  fuel  used. 


Mission  Environment  Simulation.  The  motions  and  maneuvers  of  the  craft 
are  calculated  from  the  equations  of  motion  in  a vehicle  fixed  coordinate 
system  with  origin  at  the  pilot.  However,  the  environment  external  to  the 
vehicle  is  expressed  in  an  earth  fixed  coordinate  system  with  origin  at  the 
shoreline.  These  systems  and  their  relationship  are  described  in  Section  IV. 
The  seaway  and  surf,  the  land,  and  the  ^ind  make  up  the  external  environ- 
ment inputs  to  the  vehicle.  It  is  assumed  that  the  seaway  is  made  up  of 
sinusoidal  waves  offshore  which  are  combined  with  the  vehicle-generated 
wave  system  to  model  the  total  height  of  the  water  beneath  the  craft.  The 
steepening  and  breaking  of  the  waves  as  they  approach  the  beach  is  modeled. 
Over  land  the  seaway  and  wave  calculations  are  omitted  and  ground  eleva- 
tions used  directly.  The  wind  velocity  and  direction  are  input  to  the 
simulation. 
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SECTION  11! 

MATHEMATICAL  MODEL 


INTRODUCTION 

The  mathematical  model  is  intended  for  implementation  on  a Sigma  7 
digital  computer  at  the  NAVTRAEQUIPCEN  in  Orlando,  Florida.  Digital 
to  analog,  and  analog  to  digital  interfaces  exist  for  motion,  control  and  sen- 
sor simulation.  As  the  model  has  been  developed,  it  has  been  programmed 
in  FORTRAN  IV  and  run  on  the  Draper  Laboratory  IBM  360.  It  is  assumed 
that  this  program  will  be  of  use  to  the  simulation  programmers.  It  should 
be  pointed  out,  however,  that  it  was  intended  as  a test  of  the  mathematical 
model  as  the  pieces  were  developed  and  no  attention  has  been  paid  to  effi- 
ciency or  running  time.  It  is  part  of  a mathematical  rather  than  a pro- 
gramming task  and  was  written  by  engineers  rather  than  programmers.  The 
real  programming  task,  therefore,  remains  to  be  done. 

EQUATIONS  OF  MOTION. 

BASIC  ASSUMPTIONS.  The  equations  of  motion  for  the  ACV  simulation  are 
the  standard  six  degree  of  freedom  equations  used  for  ship  motion  and  con- 
trol analysis.  They  are  fully  derived  in  the  standard  Naval  Architecture 
texts,  the  most  complete  being  Stability  and  Control  of  Ocean  Vehicles  by 
Martin  A.  Abkowitz,  MIT  Press  1969.  For  ACV  modeling  purposes  the 
origin  of  coordinates  is  taken  at  the  pilot  location  since  that  is  the  point  of 
interest  for  the  simulation.  Pitch  and  roll  angles  are  assumed  small  and 
the  appropriate  linearizations  made  in  trigonometric  functions  of  these  singles, 

The  water  and  land  surface  are  of  course  developed  in  an  earth  fixed 
level  coordinate  system  with  the  vertical  axis  truly  verticaL  This  creates 
slight  mismatches  between  heights  on  the  vehicle  and  ground  or  water  ele- 
vations. However,  for  small  roll  and  pitch  the  errors  are  negligible  and 
not  worth  the  computer  time  to  correct. 

Coordinate  Transformations.  The  two  coordinate  systems  employed  for 
the  simulation  and  the  relationship  between  them  are  described  below: 
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Definition  of  Terms 


eoj,  eo9,  eOg 
xl*  X2J  x3 


XOj,  xo2,  x©3 


xr  xj 

*°3 

x4*  X5’  Xg 


m 


II*  *2’  *3 


g 


xslc'  xs2c*  :<s3c 


FTj 


unit  vectors  forward,  to  starboard,  and 
down  on  the  craft. 

unit  vectors  North,  East  and  down. 

distances  from  ship  based  origin  (Pilot 
cabin)  in  es’1 , es"9,  es^. 

distances  from  earth  fixed  origin  (shoreline 
at  mean  sea  level)  in  eoj,  eo^* 

first  and  second  time  derivatives  of  x^  . 

time  derivative  of  xo. . 

rotations  about  esj,  es2,  es^  axes,  (roll, 

pitch,  yaw)  Xg  is  the  heading  angle  measured 

clockwise  from  North,  x.  and  x_  are  assumed 

4 o 

small. 

mass  of  vehicle. 

moments  of  inertia,  assumed  to  be  principal 
moments,  about  es"j,  eSg*  es^  axes  at  ship 
based  origin. 

acceleration  of  gravity. 

position  of  center  of  gravity  in  ship  based 
coordinates  • 

jth  component  of  total  force  vector  acting  on 
vehicle  exclusive  of  gravitational  forces  which 
are  included  in  the  equations  of  motion. 
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Coordinates 


The  transformation  from  earth  to  ship  system  is  defined  by: 


es.  = TSO-.  eo. 
i 3 

= TSOy  esj 


(1) 

(2) 


Where, 


TSOj^  = cos  xg 


TS021  = -sin  Xg 


TS031  = Xg  cosXg  +x^  sinxg 


TSOjg  = sin  Xg 


TS022  " COS  x6 


TSO32  = Xg  sinxg  - x4  cos  xg 


TSOX3  ° *5 
TS023  = x4 


tso33=  l. 


The  location  of  a point  p whose  location  on  the  ship  is  defined  by  xsj,  xs2»  xs3 
relative  to  the  pilot  cabin  is  therefore  expressed  as  follows  in  navigational 
coordinates: 


xo„,  = xo,  + xs,  cos  x_  - xs0  sin  xc  + xs0  (x_  cos  xa  + x,  sin  xe) 

pi  11  62  63o  64  0 

SOp2  = X°2  + XS1  Sin  x6  + XS2  COS  X6  + >:s3  (x5  £in  X6  " X4  CO£  X6) 

XOp3  = X°3  " XS1  X5  " XS2  X4  * XS3* 


(3) 


If  xSj  and  xs2  » xs3 


pl = 

XOj 

+ xs^ 

cos  Xg  - xs9  sin  x( 

P2  = 

xo9 

+ XSj 

sin  x_  + xs„  cos  x, 
0 Z < 

> « =5 

xo„ 

- xs. 

X_  + X50  X 4 XS_. 

p3 

3 

1 

a 2 4 3 

(4) 
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Equations  of  Motion.  The  equations  of  motion  are  integrated  using  Euler 
integration.  T5e  acceleration  of  the  pilot  location  is  represented  by: 


X.  = TFX./VMASS  - X„X  + X.X,  4 XS,  .!Xe2  4 X 
1 1 d362  1c  / 5 61 


35,  jx,X_  - X,!  - XS  »XX,  tXJ 
2c  { 4 5 6«  3c  j 4 6 5', 


(5) 


x2  = TFX2/VMASS  - x6xt  + x4x3  4 xs2c  |x62  + x42 ' 
- XS3c  !X5X6  - x4 ! -XSlc  1*5*4 


* * • • • • . - , - 

X3  = TFVVMASS  - X4X2  + X5X1  +XS3c  jX4  4 X5  j 

- XS  |x  X - X_j  - XS,  J X X_  4 X .! 
lc  f 6 4 5}  2c  i 6 5 4i 


X4  = [tFX4  - X.X6  jvi3  - VI2;  VMASS  | XS^fXj  4 X .X, 


XX)  - XS  (X  4 X.X  - X .X  } 4 XS.  x XS, 

1 . 3c  2 6 1 4 3 ic  2c 


(X4X6  - x5)  - xslc*  xs3c(X4x5  4 x6)  4 xs2c  x xs3c 
(X62  - X52  ^1 


Urr 


X = ITFX  - X.X.  )VI  - VI  I - VMASS  |XS_  (X*  4 X,  - X.X 


o 6 4 I i 3) 


I 3c 


5 3 & 2 


- XVX3  4 X,X2  - X.X,)  4 XS^*  XS3c  (X5X4  - X ^ - *Sx  XS^tX.X,  4 x4> 
4 xs3c  x XS]C(X.2  - x62)  Ij  /VI2 
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• • • • !*••••• 

X.  = TFX,  - X X | VI  - VI,  j - VM  ASS  XS,  (X,  + X .X,  - XX,) 
6 6 4 5 | 2 1)  ( lc  2 6 1 4 3 


r* 

if 


i 


- XS.  (X  + X X,  - X,X_)  + XS  x XS,  (X.X.  - XI 

2c  1 5 3 6 2 3c  lc  6 5 4 

• » ••  • * *1 

- XS  x XS.  (X,X  + X.)  + XS  x XS.  (X.  - X,  )!  I /VI 

3c  2c  6 4 5 lc  2c  5 4 j] 


The  integration  proceeds  as  follows: 


DUM,, 

= X 

N 

N 

X = 

• • * 
X,  + X 

N 

N N 

X = 

X„T  +0.5 

N 

N 

* ! 


(6) 


u 


For  the  pilot  location  in  navigational  coordinates:  (For  navigational 
purposes  heave  velocity  contributes  negligibly  to  location) 

S^  = sin  (X^)  -*■  sine  of  heading  angle 

C , = cos  (X  ) . , , , 

6 6 cosine  of  heading  angle 


i cm 


XO  ,=  X,  x C.  - X,  X S, 
1 I h c o 


X°2=  xi  * s6  + x2  X c6 


XOj  = X0X  + X0L  At 


*♦  North  pilot  velocity 
-*•  East  pilot  velocity 
New  pilot  North  of  start. 


(7) 


XO^  = XO^  + XO^At  -*■  New  pilot  East  of  start 


HPX01  = XO,  + HPXS1  x C,  - HPX32  x S, 

N 1 N 6 N 6 

New  North  location  of  position  N on  hull  for  determination  of 
seaway  height. 
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The  following  characteristics  and  their  values  used  at  present  are  below: 


IM- 


VMASS 

X1SC 

X2SC 

X3SC 

CSKRT 


ISYS 


X2R  (1,  2) 


X2N  (1,  2) 


X2P  (1,  2) 


VWIND 

AWIND 


SCAREA 
F CARE A 
LCUSH 


CHORD 

RSAREA 


mass  of  vehicle  10879.  5 slugs 
X location  of  eg  -30.0  feet 

Y location  of  eg  - 18.  0 feet 

Z location  of  eg  +8.0  feet 

skirt  stiffness  parameter  0.0112 
skirt  reaction  time  8.  0 

moment  of  inertia  of  rotating  machinery  3.  533 
X location  of  ruduers  -67. 1 feet 

Y location  of  starboard,  port  rudders  -4.  0 feet 

-32. 0 feet 

Z location  of  rudders  0. 0 feet 
X location  of  nozzles  -23. 16  feet 

Y location  of  nozzles  -0. 25  feet,  -35.  75  feet 
Z location  of  nozzles  -3.0  feet 

Y location  of  propellers  -4.  0 feet  -32.  0 feet 
Z location  of  propellers  0.0  feet 

Wind  velocity  34.  0 feet/ sec 

Wind  direction 

density  of  air  0.00237  slugs /ft ^ 

vehicle  moment  of  inertia,  about  1 axis  = 5.  672  x 10* 

r 

vehicle  moment  of  inertia  about  2 axis  = 1.629  x 10 1 

r 

vehicle  moment  of  inertia  about  3 axis  = 2.057  x 10 ' 

initial  location  of  vehicle  pilot  (North)  - feet 

initial  location  of  vehicle  pilot  (East)  - feet 

2 

planform  area  for  windage  = 3200  ft 

2 

frontal  area  for  windage  = 836  ft 
length  of  cushion  = 77  ft 
duct  area  = 123  ft^ 
duct  diameter  = 11.25  ft 
duct  chord  -•  4.  67  ft 

2 

rudder  (lifting  area  = 47.2  ft 


SDOLD  (1-— 27) 
PFAN  (1—4) 
QFAN  (—  ) 
QOC  (1—4) 
PGC  (1—4) 
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- initial  skirt  depths  - 4.  5 for  periphery,  4.0  for  dividers 

- initial  fan  pressures  - 300  psf 

- initial  fan  flows  - 7000  ft^/  sec 

- initial  flow  under  skirts  - 3000  ft  /sec  per  compartment 

- initial  cushion  pressure  - 109.  psf. 


I 

£ 

4 

§ 

% 

| 

5 

| 

1 

i 
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Forces  on  the  Vehicle.  The  force  and  moment  components  on  the 
vehicle  are  expressed  as  follows: 


TFX1 


= TPROPj  + TPROP2  - (TNOZx  + TNQZ2)  COS  (PSI)  (9) 

+ RDRAG  + XBDRAG  + XMDRAG  * WDRAGX  + SDRX  + XGRAV 


TFX2  = - (TNOZ1  + *rNOZ2)  SIN  (PSI)  + YBDRAG  + RLIFT  + YMDRAG 
+ YDUCT  + WDRAGY  + SDRY  + YGRAV 

TFX3  = FX3CT  + ZGRAV 

TFX4  = FX4CT  + ROLLA  + ROLLMD  + ROLLD  + SIN(FSI)  (TNOZj 
+ TNOZ2)  X3N  + WMX4  + RROLL  + KGRAV  + SDRR 

TFX5  = FX5CT  + PITCHA  + PITCHM  + RPITCB  - COS  (PSI) 

(TNOZj  + TNOZ2)  X3N  + (TPROPj  + TPROPg)  X3P 


TFX6 


+ WMX5  + SDRP  + MGRAV 


= - SIN  (PSI)  (TNOZj  + TNOZ2)  XIN+  COS  (PSI)  (TNOZx  x X2NX 
+ TNOZ2  x X2N2)  - X2PX  x TPROP1  - X2P2  x TPROP2 
+ YAWBD  + YAWMD  + YAWD  + RYAW  + YAWDC 


+ WMX6  + SDRYM  4-  NGRAV 


where. 


'j'FXl 

■ =Total  force  and  moment  components  in  directions  X,  to  Xc 
TFX6  1 b 

TPROPj  2=Propeller  thrusts 
TNOZf  2 =Nozzle  thrust 
PSI  = Nozzle  angle 

RDRAG  = force  on  rudder  in  eSj  , direction  (normally  negative) 
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XMDRAG 

Force  due  to  intake  air  momentum  change  in  *eSj 
direction  (normally  negative) 

XBDRAG 

= 

Aerodynamic  force  in  es*^  direction  (normally  negative) 

WDRAGX 

= 

Force  due  to  wave-cushion  pressure  interaction  in  es^ 
(normally  negative) 

SDRX 

= 

Skirt  drag  force  in  es ‘ direction  (normally  negative) 

YBDRAG 

= 

Aerodyamnic  force  in  SSg  direction 

RLIFT 

= 

Force  on  rudder  in  es2  direction 

YMDRAG 

= 

Force  due  to  intake  air  momentum  change  in  esl 
direction 

YDUCT 

= 

Force  on  ducts  in  es*2  direction 

WDRAGY 

= 

Force  due  to  wave-cushion  pressure  interaction  in  es2 
direction 

SDRY 

= 

Skirt  drag  force  in  es2  direction 

FX3CT 

= 

Force  on  vehicle  due  to  cushion  pressure  in  es^  direction 

FX4CT 

= 

Moment  on  vehicle  due  to  cushion  pressure  in  roll 
direction 

ROLLA 

= 

Roll  moment  due  to  aerodynamic  forces 

ROLLMD 

= 

Roll  moment  due  to  intake  air  momentum  change 

ROLLD 

= 

Roll  moment  due  to  duct  side  force 

X3N 

= 

Vertical  position  of  nozzle 

WMX4 

a 

Roll  moment  due  to  wave-cushion  pressure  interaction 

♦ 


23 


NAVTRAEQUIPC12N  73-C-0138 


SDRR  = Roll  moment  due  to  skirt  drag 

FX5CT  = Pitch  moment  due  to  cushion  pressure 
* 

PITCHA  = Pitch  moment  due  to  aerodynamic  forces 
PITCHM  = Pitch  moment  due  to  intake  air  momentum  change 
R PITCH  = Pitch  moment  due  to  forces  on  rudder 
X3P  Longitudinal  position  of  propellers 

WMX5  = Pitch  moment  due  to  wave-cushion  pressure  interaction 
SDRP  = Pitch  moment  due  to  skirt  drag 
YAWBD  = Yaw  moment  due  to  aerodynamic  forces 
YAWMD  = Yaw  moment  due  to  intake  air  momentum  change 
YAWD  = Yaw  moment  due  to  duct  sideforce 


RYAW  = Yaw  moment  exerted  by  rudders 


YAWDC  = Yaw  damping  coefficient 


WMX6  = Yaw  moment  due  to  wave-cushion  pressure  interaction 


SDRYW  = Yaw  moment  due  to  skirt  drag 


RROLL  = Roll  moment  due  to  rudder  forces 


XGRAV  = Gravity  force  in  es  direction  (-mgX  ) 

X D 

YGRAV  = Gravity  force  in  es  direction  (mgX  ) 

ZGRAV  = Gravity  force  in  es^  direction  (mg) 

KGRAV  = Gravity  moment  in  es^  (ZGRAVxX2SC-YGRAVxX3SC) 


MGRAV  = 


Gravity  moment  in  es&  (XGRAVxS3SC-ZGRAVxXlSC) 


Gravity  moment  in  es  (YGRAVxXlSC-XGRAVxXlSC). 
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MODELING  OF  EFFECTORS  AND  CONTROLS 

EFFECTOR  CONTROL  LOGIC.  Angle  Control  Logic.  For  controls  acti- 
vated by  the  pilot  the  routine  STEER  is  used  to  compute  the  response  to 
pilot  input.  STEER  works  as  follows: 


INPUT 
C.  A.  M 

T — 

DUM  = 1.1  X A X At 

— ..nrzz: 

YES  iS  NO 

1^  - Cl  < DUM  I 


£ = PRESENT  POSITION  OF  DEVICE 
C = ORDERED  POSITION  OF  DEVICE 
A = VELOCITY  OF  DEVICE 
At  = TIME  INTERVAL 


£ = C 


RETURN 


A|^  = 

A At 

IFC<* 

A^  = — A^/ 

^ 

-a. 

<1 

+ 

RETURN 

The  actuator  ( 4*  ) moves  toward  its  ordered  70c  it  ion  (C)  at  rate  (A). 

The  nozzle  angle  PSI  is  controlled  in  two  modes.  If  the  mode  SWITCH  is 
positive,  angles  between  -90  and  +90  are  intended,  and  the  thrust  is  aft. 


SWITCH 


BOW 

♦ k*  / 


STARBOARD 


SWITCH  - 


-THRUST 


THRUST 
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If  the  mode  switch  is  negative,  angles  between  90°  and  270°  are  intended,  and 
the  thrust  is  forward.  The  ordered  control,  besides  switch,  is  the  steering 
wheel,  which  operates  between  -30  and  +30  . The  steering  wheel  angle  is 
ASTEER,  When  SWITCH  is  positive  3 x ASTEER  is  the  ordered  angle.  When 
SWITCH  is  negative,  3 x ASTEER  + 180°  is  the  ordered  ang1e.  The  logic  is 
handled  by  a dummy  variable  PS  which  represents  the  nozzle  angle  PSI  in  the 
SWITCH  + mode,  and  the  nozzle  angle  less  180°  in  the  SWITCH  - mode. 
Therefore,  the  logic  is: 

PS  = PSI  (11) 

If  SWITCH  is  negative,  PS  = PSI  -180° 

call  STEER  (PS,  3 x ASTEER,  ACONST,  DELT) 

If  SWITCH  is  negative,  PSI  = PS  + 180° 

ACONST,  the  angular  velocity  of  the  nozzle,  is  50°/  sec. 


For  the  rudder  angle  RANGLE,  the  input  from  the  pilot  is  RCONT, 
the  desired  rudder  angle  from  the  foot  pedal  controls.  STEER  is  caUed  with: 

(RANGLE,  RCONT,  RCONST,  DELT).  RCONST,  the 
rudder  rate,  is  30°/ sec. 

For  the  propeller  pitch  angles,  the  input  from  the  pilot  is  APROP(l) 
and  (2),  the  starboard  or  port  pitch  lever  positions  STEER  is  called  with: 

(BETA  (N),  APROP  (N),  PCONST,  DELT).  PCONST,  the 
pitch  rate,  is  20  /sec. 

Machinery  Control.  Variable  definitions:* 


NIDES(l) N1DES(6)- 


N2GOV(l) N2GOV(2)- 


Desired  gas  producer  control  rpm 
set  by  pilot  on  each  turbine. 

Maximum  power  turbine  control  rpm 
set  by  pilot,  starboard  and  port. 


Nl(l) Nl(6) 


Actual  gas  producer  rpms  on  each 
turbine. 


N2(l) N2(2) 


SlCON 


Actual  power  turbine  shaft  rpm,  star- 
board and  port. 

Time  c<pstant  of  N1  control  (input, 

. 3 sec”1). 


S2CON- 


Time  constant^  N2  governor  control 
(input,  .5  sec  /. 


NAVTRAEQUIPCEN  73-00138-1 


Machinery  Control  Logic.  The  pilot  controls  on  the  turbines  consist  of  a 
desired  fuel  flow  rate  on  each  turbine  (N1DES)  and  a governing  (upper  limit) 
control  on  the  starboard  and  port  power  turbine  output  rpm  (N2GOV).  If 
N2GOV  is  exceeded  the  corresponding  fuel  rates  are  decreased. 

Logic  J = 1 (for  starboard  side),  n = 1, 2, 3 

J * 2 (for  port  side),  n = 4,  5, 6 

if  N2..  < N2GOVj  (12) 

N1  = N1  + Si  CON  jNlDES  -N1  t At 
n n ' n n' 


if  N2_.  >N2GOVj 

N1  = Nl  + S2CON  1n2GOV  - N2  (At. 
n n 1 n n 

Propeller  Modeling,  The  thrust  and  torque  characteristics  of  the  prop^ers 
are  based  entirely  on  the  Bell  horizontal  plane  simulation.  This  is  described 
in  their  Preliminary  Design  Summary  Report,  appendix  to  paragraph  2. 2. 2, 
"Three-Degree-Of-Freedom  Maneuverability  and  Control  Simulation. 11  Many 
of  the  effector  and  control  descriptions  in  the  present  mathematical  model  are 
taken  from  this  Bell  Aerospace  simulation.  In  the  notation  of  the  present 
simulation: 


For  propeller  thrust: 


if  0j  > 0. 0 


TPROPj  = |338.  Pi  + 4.36  0j2-O.  1715  XWIND2  - 1.43  ^jXWIND 
|NPROPj/l250|  2 (13) 


if  pi  <0.0 

TPROPj 

For  propeller  power: 


SO. Pi  - . 1715  X WIND2J  J NPROPj/ 1250 j 2. 


if  0j  >10.,  Cl  =0.0 

if  Pi  < 10.,  Cl  = 5o|xwIND/127.5|  jio  - 03  | 

HPROPj  = (NPROPj/ 1250.  )3  1 450.  + 23.05  0j  + 2. 56  0 j2  +C1 


where 

I3j 

3 

TPROP. 

XWIND3 


Pitch  angle  of  propeller  j 
1 for  starboard,  2 for  port 
Thrust  in  pounds  on  propeller  j 
Apparent  head  wind  velocity 


(14) 


NPRCP. 

3 

HPROP. 

3 
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= RPM  of  propeller  j 

= Horsepower  absorbed  by  propeller  j. 


The  forces  and  moments  exerted  by  the  propellers  are: 


in  ESj  direction 


in  yaw  direction 


2 

£ 

3=1 

2 

L 

3=1 


TPROP. 

3 


- TPROP.  x X2P.. 
3 3 


(15) 


Thrust  Nozzles.  The  thrust  exerted  by  ?acb  nozzle  is  derived  directly  from 
the  momentum  flux  of  the  air  exciting  from  the  nozzle.  The  effect  of  rela- 
tive wind  velocity  on  this  force  is  considered  small. 


TNOZ.  = 2.44x  10"4QNOZ2 
3 3 


(16) 


i-  \ 

W 


where 


TNOZ  = 
3 

QNOZ_.  = 


Thrust  in  pounds  on  nozzle  j 
1 for  starboard,  2 for  port 
Flow  through  nozzle  j (ft/ sec). 


The  flow  through  the  nozzle  is  obtained  directly  from  the  fan  pressure: 


O 


QNOZ.  = -346.  '/jPFAN.1  sgn  (PFAN.)  (17) 

3 3 3 

where 

PFAN . is  the  pressure  developed  by  fan  j in  psf . 

3 

Rudders.  The  forces  exerted  on  the  rudders  are  based  on  the  Bell  PDSR 
simulation.  First,  the  velocity  at  the  rudders  including  vehicle  motion, 
effect  of  the  ducted  propellers,  and  wind  must  be  found.  The  square  of  the 
axial  velocity  at  rudder  j is: 

If  TPROP.  > 0 
3 

VDUCT2.  = (VA  cos £aa)2  + TPROP./ (PAE)  (18) 

3 3 - 

If  TPROPj  < 0 

VDUCT2.  = (VA  cos^aa)2  + TPROP./(2PAE) 

• 3 
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where, 

VA 

)3aa 

TPROP. 
. 3 3 

AE 


Apparent  wind  velocity  in  ft/  sec 
Apparent  wind  angle  in  radians 
Thrust  on  propeller  j 
1 for  starboard,  2 for  Dort 
Air  density  in  ^Lugs/ft*^ 

Duct  area  in  ft  (123  ft  ). 


The  dynamic  pressure  available  at  the  rudders  is  then 


PDUCT. 

3 


1/2  PVDUCT2.. 
' 3 


(19) 


The  rudder  lift  and  drag  coefficients  for  a rudder  angle,  RANGLE  (in  degrees) 
are: 


CLEFT  = . 053  RANGLE  /0. , 

-3  2 

CDRAG  = . 02  + .422  x 10  RANGLE 

in  the  small  angle  of  attack  range.  Stagnation  is  assumed  to  occur  at  20°: 

If  RANGLE  >20°,  CLIFT  =1.06. 

The  longitudinal  and  side  forces  due  to  the  rudders  are  then: 

RDRAG  = -CDRAG  x P.SAREA  x (PDUCTj+PDUC^)  (21) 

RLIFT  = CLIFT  x RSAREA  x (PDUCT^PDUCT^ 

2 

where  RSAREA  is  the  rudder  area  (47. 2ft). 

Due  to  the  location  of  the  rudders  relative  to  the  pilot  location  origin  yaw, 
roll  and  pitch  moments  are  developed: 

RYAW  = CDRAG  x RSAREA  (X2R  xPDUCT  +X2R_ 

1 1 ^ 

x PDUCT2)  -f  XI R x RLIFT 
RROLL  = -X3R  x RLIFT 

RPITCH  = +X3R  x RDRAG. 
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Engine  Modeling.  Basis  of  Model.  The  turbine  characteristics  are  based  on 
the  Bell  Aerospace  data  shown  in  figure  6.  A curve  fit  is  used  to  find  the 
optimum  power  turbine  rpm,  Ng,  for  gas  generator  rpm  N^.  The  gas  gene- 
rator rpm  N has  been  set  by  turbine  control,  as  shown  above.  The  shaft 
horsepower  generated  at  optimum  Ng  is  found  from  a curve  fit  to  figure  6 and 
the  actual  shaft  horsepower  from  a parabolic  fit  to  the  off-optimum  loss  for 
the  actual  N„.  The  torques  required  by  fans  and  propellers  are  then  com- 
puted at  turmne  rpm  Ng  and  the  torque  put  out  by  the  turbin--  •»  compared 
with  that  absorbed.  Any  excess  accelerates  the  system. 

Variable  definitions. 

Nl(l) Nl{6)  —Gas  generator  rpm  on  each  turbine.  Turbines 

1 to  3 are  starboard. 


N20PT(1) N20PT(6) 

SHPOPT(l) SHPOPT(6) 

SHP(l) SHP(6) 

GPH(l) GPH(6) 

ETOR(l) ETOR(6) 

TETOR(l) TETOR(2) 

N2(l) N2(2) 

n,  j 


NFAN(l) N.FAN(2) 

NPROP(l) NPROP(2) 

HPFAN(l)  — HPFAN(2) 

PFAN(l) PFAN{2) 

QFAN(l) QFAN(2) 


—Optimum  power  turbine  rpm  N2  for  each 
turbine  Nl. 

—Horsepower  generated  by  each  turbine  if 
operated  at  optimum  N2. 

—Actual  horsepower  generated  by  each  turbine. 

—Gallons  per  hour  of  fuel  used  by  each  turbine. 

—Torque  produced  by  each  turbine  in  foot  pounds. 

—Total  turbine  torque,  starboard  and  port. 

—Power  turbine  rpm.  N2(l)  is  starboard,  N2(2) 
is  port. 

—Dummy  index  variables.  For  engines  n - 1.2,3 
are  starboard  and  correspond  to  j=l.  n=4,5,6 
are  port  and  correspond  to  j=2. 

—Starboard  and  port  cushion  fan  rpm. 

—Starboard  and  port  propeller  shaft  rpm. 

—Starboard  and  port  power  absorbed  by  cushion 
fans. 

—Starboard  and  port  fan  pressures  from  cushion 
air  dynamics. 

—Starboard  and  port  fan  airflows  from  cushion 
air  dynamics. 
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NETAF 

TNOZ(l) 

BETA(l) 

XWIND 


— Fan  efficiency. 

TNOZ(2)  —Thrust  from  nozzles. 

BETA(2)  — p,  propeller  pitch  angle  in  degrees. 

—Relative  head  wind  velocity  in  direction 
from  Main  program. 


TPROP(l) TPROP(2)  —Starboard  and  port  propeller  thrust. 

HPPROP(l) HPPROP{2)— Horsepower  absorbed  by  starboard  and  port 

propellers. 

FTOR(l) FTOR(2)  —Torque  used  by  fans  at  N2  rpm,  -starboard 

and  port. 

PTOR(l) PTOR(2)  —Torque  used  by  propellers  at  N2  rpm,  star- 

board and  port. 

DN2DT(1)  — DNDT(2)  — Turbine  power  shaft  accelerations  in  radians/ 

sec  . 

ISYS  —Equivalent  moment  of  inertia  of  machinery 

as  soon  from  a turbine  power  shaft. 

Machinery  Characteristics.  The  optimum  turbine  rpm  N2  is  determined 

for  each  turbine  for  its  gas  generator  rpm  Nl  from:  (Curve  fit  from  figure  4). 

(22) 


N20PT.  * 11672  - 1.S0457  Nl  - + 1.21488  x 10~4  Nl2. 
n n n 


The  shaft  horsepower  each  turbine  would  generate  at  this  optimum  N2  is 
then  determined: 

SHPOPT  = 14684.  - 2.3976  Nl  + 9.  796  x 10"5  Nl2. 
a n n 

The  actual  horsepower  produced  by  each  turbine  is  then  calculated: 


(23) 


SHP  = SHPOPT  (o.l  * 1.8  (N2j/N20PT-)  -0. 9(N2j/N20PT« (24) 
nl  n n J 

The  fuel  consumption  rate  is  calculated  for  each  turbine  from  the  gas  gene- 
rator rpm 


GPH  = 757.  6377  - . 121784  Nl  + 5.  20379  x 10~6N12. 
n n n 


u 


(25) 
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The  torque  is  then  determined  for  each  turbine: 

ETORn  = 3300  SHPn/(27rN2J). 

The  total  starboard  torque  is : 

TETOR(l)  » ? E-TOR 
1-3  n 

and  port: 

TETOR<2)  = * ETOR  . 

4-6  n 

The  gear  ratios  in  the  shafting  system  are  then  used  to  find  propeller  and 
fan  speeds: 


NFANj  = .1297  N2; 


3 = 1>2 


NPROPj=  . 6427  NFANj  j = 1,  2. 

The  fan  power  absorbed  is  then  calculated  from  the  efficiency  v , 

HPFANj  = PFANj  QFANj/  (540  NETAF). 

From  the  fan  and  propeller  powers  the  torque  absorbed  from  the  turbines 
derived: 

FTORj  = 33000  HPFANj  / (2irN2j) 


PTORj  = 


33000  HPROPj  / (27rN2j). 


The  difference  between  the  torque  produced  by  the  turbines  (TETOR(l)  and 
TETOR(2)  ) and  that  absorbed  by  fans  and  propellers  loads  to  acceleration 
of  the  shafting: 

DN2DTj  = {teTORJ  - FTORj  - PTORj  } / ISYS.  (5 

The  shaft,  rpms  are  then  updated: 

N2j  = N2j  + DN2DTj  At.  (3 
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CUSHION  SYSTEM  MODEL 

The  determination  of  cushion  pressure  is  governed  by  the  fan  pressure- 
flow  relationship,  the  change  of  cushion  volume  with  time  due  to  waves  or 
craft  motion,  and  the  pressure-flow  relations  beneath  the  skirts. 

For  speed  in  computation,  internal  functions  are  set  up  for  computation 
of  individual  cushion  volumes  and  skirt  gap  areas  depending  on  the  vehicles' 
present  relation  to  the  water  or  ground  surface.  These  routines  are  described 
below.  At  each  time  step  the  volume  of  each  air  compartment  is  determined 
using  the  internal  function  HTAV  and  the  height  of  the  bottom  plating  above 
water,  HOW,  at  the  35  points  shown  in  figure  7.  The  depths  of  the  skirts 
below  this  bottom  plating  is  a function  of  cushion  pressure  and  their  determi- 
nation is  described  below.  The  gap  between  skirts  and  water  is  calculated 
using  the  heights  above  water  HOW  and  skirt  depths  SD  at  the  peripheral 
points.  These  gaps  are  used  with  the  internal  functions  QL  and  QS  to  deter- 
mine skirt  pressure-flow  relationships.  This  volume  and  gap  characteristic 
determination  is  described  below. 

With  geometry  determined,  the  pres  sure -flow  relations  are  solved  and 
the  re  malting  forces  on  the  vehicles  are  derived  below. 

Integration  formulas.  For  the  volume  under  each  of  the  four  cushion  com- 
partments an  integration  formula  is  established  for  the  average  height  under 
the  following  configuration:  (One  of  four  compartments) 


1( 

y 


1 1 

11 

f-12  D 

1 i 

t 

i 

3 8 1 

2 

Assuming  cubic  curves  in  x and  parabolic  in  y the  internal  function  HTAV  is 
established  for  the  height  for  volume  calculations: 

HTAV  (AN1,  AN2,  AN3, AN12)  = .0203333  {AN1  + AN2  + AN 3 + AN4}  (33) 

+ .0625  {AN5  + AN6  + AN7  + AN8}  + .08333  {AN9  + AN10  + .25  AN11  + AN12}, 

The  height  of  the  hull  bottom  plating  above  the  water  or  ground  surface 
are  calculated  for  each  cushion  aad  the  result  multiplied  by  the  cushion 
planform  area  to  obtain  compartment  air  volume  as  a function  of  vehicle 
orientation  and  water  or  ground  surface,  ‘»?(x  , y , t). 
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For  calculation  of  pressure  differences  across  the  skirt  gap  a discharge 
coefficient  of  0. 6 is  assumed: 


AP  = 


A2  2 CT 


where  Q is  the  flow  rate  under  the  skirts  and  A is  the  gap  area  under  the 
skirt  for  this  compartment. 

Assuming  skirt  perimeter  configurations  for  each  compartment  as  follows: 


functional  routines  are  established  for: 


A v 2cd/p 


In  X a cubic  gap  is  assumed: 


QLtANSl. . . ANS4)  = 94. 16  {ANSI  + 3xANS2  + 3xANS3  + ANS4}  . 

In  Y a parabolic  gap  is  assumed:  (36) 

QS(ANS1,  ANS2,  ANS3)  = 62.  76  {ANSI  + 4xANS2  + ANS3{  . 

As  the  heights  above  water  or  ground  of  the  appropriate  points  on  the 
cushion  shift  periphery  are  input,  the  quantity 


QL  or  QS 


. /2C 
-4 Vp-^ 


results. 


Determination  of  Cushion  Volumes  and  Skirt  Gaps. 


Cushion  Volumes.  The  heights  of  each  point  on  the  bottom  plating  are 
determined  measured  from  mean  sea  level: 

HHULLn  = -X3  -f  X5  (HPXSln)  -X4(HPXS2n)  - HPXS3n).  (37) 


m 


HOW 


« HULL  - ETA  . 
n n 


where 
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HPXSlr 

HPXS2r 

HPXS3. 


Is  the  location  of  point  n as  measured  from 
the  pilot  in  es^,  es2,  es^  . 


ETA_ 


Is  the  height  of  water  or  ground  at  point  n 
above  mean  sea  level. 


The  volume  of  each  cushion  compartment  (VOLC  ) is  then  determined  with 
the  heights  above  and  the  integration  formula  HTAV. 

The  change  of  cushion  volume  with  time  may  then  be  calculated.  This 
is  usually  called  "wave  pumping"  and  is  the  flow  into  the  cushion  required 
by  its  rate  of  change  of  volume. 

QPUMP.  = (VOLC.  - VOLC.  (old)) /At.  (39) 

3 3 3 

Skirt  Gaps.  An  inner  loop  is  set  up  for  the  calculation  of  the  depth  of  the 
finger  skirt  periphery  below  the  hard  hull.  The  skirts  are  assumed  to 
respond  to  pressure  variations  in  the  compartment  according  to  the  following 
logic: 

Y = 109.  - PGC  difference  between  equilibrium  value  and 

n n present  compartment  air  pressure. 


> 66.  9 Y_  = 66.  9 


limits  on  skirt  motion 


Y <-66. 9 Y = -66. 9. 
n n 

SDPRO  = 4.  5 + CSKRT  xY  - 8.  325  x 10-7  Y3  (^ 

n n n 

cubic  fit  to  steady  skirt  depth  as  a function  of  pressure. 


SDDT  = (SDPRO  - SDOLD  ) x TC 
n n n 

approach  to  equilibrium. 


skirt  depth 


SDDT  = SDOLD  + SDDT  x DELTAT  ^new  skirt  depth.  (43) 
n n n r 

The  height  of  the  skirt  perimeter  above  (or  below)  water  is  calculated, 
SAWn  = HOWn  - SDn.  (44) 

If  SAW  is  negative,  the  skirt  is  touching  water  at  point  n and  SAW  will 
be  used  for  skirt  drag.  Otherwise  CLR  , the  airgap  at  this  point n 


E5 
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is  calculated: 


CLRn  = SAW  (45) 

If  CLR  < 0,  CLR  = 10"6 
n n 

The  clearance  (CLR)  is  then  used  with  functions  QL  and  QS  to  determine  the 
pressure-flow  relation  below  the  skirts. 

Cushion  Pressure  Determination.  The  mathematical  model  developed  for 
pressures  ana  flow  rates  through  the  air  cushion  vehicle  support  system  is 
based  on  assumptions  of  incompressible  fluid  flow  and  on  discharge 
coefficients  for  flows  through  ducts  and  openings.  In  addition,  inertial 
characteristics  of  the  airmass  are  assumed  to  be  small. 

These  assumptions  greatly  reduce  the  complexity  of  the  mathematical 
model,  but  solution  of  the  resulting  state  equations  still  entails  the  solution 
of  six  simultaneous  nonlinear  equations.  These  six  equations  are  airflow 
continuity  equations  for  the  six  ''nodes"  in  the  flow  system  as  shown  in  figure  8. 
The  sum  of  all  the  flows  into  each  of  the  two  "manifolds"  downstream  of  the 
fans  must  equal  zero.  The  flow  rates  into  the  manifolds  depend  on  the  fan 
rpm  and  the  pressure  difference  between  manifold  and  atmosphere.  Like- 
wise flow  rates  out  are  related  to  pressure  differences  between  manifold 
and  cushions. 

Similar  flow  continity  equations  are  written  for  each  of  the  four  cushions, 
where  attention  must  be  given  to  the  rate  of  change  of  cushion  volume  due  to 
vehicle  motions  and  motion  of  the  water  surface.  In  figure  8 these  effects 
are  represented  symbolically  as  a piston  in  each  cushion. 

For  purposes  of  matrix  manipulation  in  this  section  all  pressures  are 
stored  in  vector  P: 

Pj  Pressure  in  cushion  1 = PGC  (1) 

Pg  Pressure  in  cushion  2 = PGC  (2) 

Pg  Pressure  in  cushion  3 = PGC  (3) 

P4  Pressure  in  cushion  4 = PGC  (4) 

P5  Starboard  fan  pressure  = PFAN  (1) 

PR  Port  fan  pressure  = PFAN  (2). 


Defining: 


QlNC(i) 
QNOZ(i) 
QFAN  (i) 
QPUMP(i) 
QlC(i) 

Q(i) 


flow  into  cushion  i from  manifold 

flow  into  manifold  i through  thrust  nozzle 

flow  into  manifold  i through  fan 

change  of  cushion  volume  i with  time 

flow  into  cushion  i from  cushion  (i-1) 

flow  into  cushion  i under  its  peripheral  skirt. 
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Cushion  Pressure  Schematic 
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The  six  continuity  equations  are: 

For  the  manifolds: 

-QINC(l)  - QINC(2)  + QNOZ(l)  + QFANU)  = 0 

(46) 

-QINC(3)  - QINC(4)  + QNOZ(2)  + QFAN(2)  = 0. 

For  cushions: 

QINC(l)  + QIC(l)  - QIC(2)  + QPUMP(l)  + Q(l)  = 0 (47) 

QINC(2)  + QIC(2)  - QIC(3)  + QPUMP(2)  + Q(2)  = 0 

QINC(3)  +QIC(3)  - QIC(4)  + QPUMP{3)  + Q(3)  = 0 

QINC(4)  + QIC(4)  - QIC(l)  + QPUMP(4)  + Q(4)  = 0. 

Most  of  these  flows  may  be  expressed  in  terms  of  the  pressures  at  the  6 
nodes  and  the  cross  sectional  area  of  the  passages  connecting  them.  The 
only  exceptions  are  the  fan  flows  which  also  depend  on  the  fan  rpms.  The 
general  form  is: 

2 — 

Q-  (area  sgn  (Ap).  (48) 

The  resulting  expressions  are: 


QINC(1) 

= 

589.  n4p(5)  - P(l)( 

sgn  (P(5)  - P(l)) 

QINC(2) 

- 

583.  >(5)  - P(2)| 

sgn  (P(5)  - P(2)j 

QINC(3) 

= 

589.  v(P(6)  - P(3)! 

sgn  (P(6)  - P(3)) 

QINC(4) 

= 

589.  v<P(6)  - P(4)l 

sgn  (P(6)  - P{4)) 

QNOZ(l) 

= 

-346.v<P(5)i 

sgn  (P(5)) 

QNOZ(2) 

= 

-346,v<P(6l. 

sgn  (P(6)) 

QIC(l) 

= 

675.  VlP(4)  - P(l)l 

sgn  <P(4)  - P(l)) 

QIC(2) 

= 

338.  V<P(1)  - P(2)| 

sgn  (P(l)  - P(2)) 

QIC(3) 

= 

675.  ,/tP(2)  - P(3)l 

sgn  (P(2)  - P(3)} 

QIC(4) 

= 

338.  v1P(3)  - P(4)l 

sgn  (P{3)  - P(4)) 

l J 
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Q(l)  = -area  (1)' (14. 5)  yiPO)T  sgn  (P(l)) 

Q(2)  = -area  (2)-  (14.  5)  v/|P(2)i  sgn  (P(2)) 

Q(3)  = -area  (3)-<14.  5)  /jp<3fc  sgn  (P(3)) 

Q(4)  = -area  (4)- (14.  5)  ,/lP(4)i  sgn  (P(4)). 

The  fan  characteristics  are  represented  by  the  following  expressions: 

QFAN  (1)  [-1280.  v/iP(5)  ••  300 J sgn  (P5)  - 300. ) (49) 

-31.  6(P(5)  - 300.)]  NF AN(1) / 2000. 

QFAN<2)  [-1280.  /|P(6)  - 300.1  sgn  (P{6)  - 300. ) 

-31. 6(P(5)  - 300.)]  NFAN(2),  2000. 

The  four  QPUMP  values  are  determined  by  the  motion  of  the  craft 
and  water  surface  and;  thus,  form  an  inhomogeneous  part  of  each  cushion 
equations. 

Substitution  of  the  above  into  equations  gives  six  simultaneous 
nonlinear  algebraic  equations  for  P(i)  given  QPUMP,  NFAN,  and  the  escape 
area  beneath  the  skirts.  These  equations  are  solved  iteratively  by  a multi- 
dimensional Newton -Raphson  procedure. 

Symbolically  the  equations  may  be  written: 
lP,  1 W.l 


IM  . ’*| 
Ik)  i ) 


where  [O^j]  represents  the  six  nonlinear  operators  on  the  vector  P and  W 
is  a vector  containing  the  inhomogeneous  terms.  If  we  introduce  an  error 
vector  E,  which  should  be  zero  when  the  equations  are  satisfied,  we  have: 


El) 

(pl) 

(wi) 

: ( 

=[°nJ  : 

•r 

E ' 
0 

(V 

(w6) 

«■'«  « W.  il  i IV I,  iWwH! 
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With  approximate  values  for  the  vector  P,  a linearization  is  possible  giving: 


iEl] 


iEl)( 

III 


pdE-l  (Pl“Plguess) 

+ Li  : (.  (52) 

LdPiJ  (p„ -p, 1 


Eg  at  P = Pguess 


6guess 


by  setting  E = 0 an  updated  approximation  to  P is  obtained  by  solving  the 
above  set  of  coupled  linear  algebraic  equations.  This  is  iterated  untL. 
suitable  accuracy  is  obtained.  In  general  one  or  iwo  iterations  are  required. 

The  derivation  of  the  error  vector  E and  the  matrix  dE^ dP.  is  as  follows: 

E(l)  = QPUMP(l)  + 5.89x102ASQRT(P(5)  - P(l)  ) + 6.  75xl02ASQRT(P(4)  - P(l)  ) 

2 (5 

- 3.38x10  ASQHT(P(1)  - P(2)  ) - 14.5  AREA(l)  ASQRT(P(1>  ) 


J ! 


dE(l)/dP(l) 


-5.89xl02  DASQRT  (P(5)-P{1)  ) - 6.75xl02  DASQRT  (P(4)-P(l)  ) 
-3.38xl02  DASQRT  (P(l)-P(2)  ) + AREA(l)  DASQRT  (P(l)  ) 


dE(l)/dP(2)  = 3.38x10  DASQRT  (P(l)-P(2)  ) 


dE(l)/dP(4) 


6.  75x10  DASQRT  (P(4)-P(l)  ) 


dE(l)/dP(5)  = 5.89x10^  DASQRT  (P(5)-P(l)  ) 


dE(l)/dP(l)  = -dE(l)/dP(5)  - dE(l)/dP(4)  - dE(l)/dP(2) 

-14.5  AREA(l)  DASQRT  (P(l)  ) 


QPUMP<2)  + 5.  89xl02  ASQRT  (P(5)-P(2)  ) + 3.38xl02  ASQRT  (P{1) 
-P(2)  ) - 6.75xl02  ASQRT  (P(2)-F(3)  ) - 14.5  AREA(2)  • 
ASQRT  (P{2)  ) 


dE{2)/dP(l)  = dE(l)/dP(2) 

dE(2)/dP(3)  = 6. 75xl02  DASQRT  (P(2)-P(3)  ) 

dE(2)/dP<5)  = 5.89xl02  DASQRT  (P{5)-P<2)  ) 

dE(2)/dP(2)  = -dE(2)/dP(5)  - dE(2)/dP(l)  - dE(2)/dP(3) 


-14.5  AREA(2)  DASQRT  <P{2)  ) 
42 
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where: 

ASQRT  (X) 

DASQRT  (X) 

E(3)  = QPUMP(3)  + 5.  89xl02  ASQRT  :<P(6)-P(3)  ) + 6. 75xl02  ASQRT  (P(2)  - 

-P(3)  ) -3.38xl02  ASQRT  (P(3)-P(4)  ) -14.5  AREA  (3)  ASQRT  <P(3)  ) 

dE(3)/dP(2)  = dE(2)/dP<3) 

dE(3)/dP(4)  = 3.38xl02  DASQRT  (P(3)-P(4)  ) 

dE(3)/dP(6)  = 5. 89xl02  DASQRT  (P(6)-P(3)  ) 

dE(3)/dP(3)  - -dE(3)/dP(6)  - dE(3)/dP(2)  - dE(3)/dP(4) 

-14.5  AREA«3)  DASQRT  (P<3)  ) 

E{4)  = QPUMP(4)  + 5. 89xl02  ASQRT  (P(6)-P(4)  ) + 3.38xl02  ASQRT  (P(3) 

-P(4)  ) - 6.75xl02  ASQRT  (P(4)-P(l)  ) - 14.5  AREA{4)  ASQRT  (P(4)  ) 


4tKi  sign  (X) 

1 = £4x 

2>|xi  dx 


dE(4)/dP(l) 

dE(4)/dP(3) 

dE(4)/dP(6) 

dE(4)/dP(4) 


E(5) 


dE(5)/dP(l) 

dE(5)/dP(2) 


dE(l)/dP(4)  • 
dE(3)/dP(4) 

5.89xl02  DASQRT  (P(6)-P(4)  ) 
dE(4)/dP(l)  - dE(4)/dP(3)  - dE(4)/dP(6) 

-14.5  AREA  (4)  DASQRT  (P(4)  ) 

-5.89xl02  ASQRT  <P(5)-P(1)  - 5. 89xl02ASQRT  (P(5)-P(2>  ) 

-3.46xl02ASQRT  (P(5)  ) -(NFAN(l)/2000.  )2  Jl.28xl03ASQRT 
<P(5)-300)  + 3.16x10*  (P(5)-300)j 
dE(l)/dP(5) 
dE(2)/dP(5) 
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dE(5)/dP(5)  = -dE(5)/dP(l)-dE(5)/dP(2)  - 3.46  x 102  DASQRT  (P5)) 

-(NFAN(l)/2000)2  Jl.  28  x 103  DASQRT  P(5)  - 300)! 

E(6)  = -5.  89  x 102  A SQRT  (P(6)  - (P(3))  - 5.  89  x 102ASQRT  (P(6)-P(4) 

-3.46  x 102  A SQRT  <P6»  - (NFAN(2)/2000.  )2  jl.  28  x 103 
A SQRT  (P(6)  -300)  + 3. 16  x 10'  (P6)  - 300)]  . 

Forces  and  Moments  due  to  wave  system.  Determination  of  the  hydrodynamic 
forces  acting  on  the  craft  is  made  by  integrating  the  cushion  pressure  times 
the  unit  vector  normal  to  the  water  surface  below  the  craft.  Since  this  surface 
tends  to  slope  down  from  bow  to  stern,  there  will  generally  be  a force  opposing 
the  craft’s  forward  motion.  Calculation  of  the  overall  shape  of  the  water 
surface  involves  both  wind  generated  ocean  waves  and  shipwaves  resulting 
from  the  craft’s  motion.  These  are  treated  in  a later  section.  Under  the 
assumption  of  uniform  pressure  within  each  chamber,  this  integral  simplifies 
in  evaluating  the  differences  in  water  surface  heights  at  opposite  ends  of 
chambers. 

The  hydrodynamic  forces  can  then  be  evaluated  as  follows  with  ETA(i) 
referring  to  the  ith  planform  point.  ETA.  is  the  sum  of  the  water  amplitude 
at  point  i due  to  the  seaway  and  the  vehicle  generated  waves: 

ETAj  = ETA.  (seaway)  + ETA.  (motion)  (54) 

Wave  Drag  and  Sway  force,  and  Yaw  Moment  from  Elevations. 

In  the  Ahead  Direction  (XI): 

W.m  = -PGC(l)  * 20  * (ETA(l)  + ETA(2)  + ETA(3)  - ETA(6)  - ETA{7)  (55) 

- ETA(8))/3 

IV D1  = WD1  - PGC(2)  * 20  * (ETA(6)  +ETA(7)  +ETA(8)  - £TA{21)  - ETA(27) 

- ETA(26))/3 

WD1  = WD1  - PGC(3)  * 20  * (ETA  (8)  + ETA(l6)  + ETA{15)  - ETA(iv) 

- ETA  (20)  - ETa(21))/3 

WDl  = WD1  - PGC(4)  * 20  * (ETA(l)  4 ETA  (11)  4 ETA(12)-  ETA(8) 

- ETA(16)  - ETA  (15))/3 


In  Ihe  Side  Direction  (X2): 

WD2  = -PGC(l)  * 40  * (ETA  (3)  4 ETA(4)  4 ETA(5)  4 ETA(6)  - ETA  (1) 
-ETA(10)  - ETA(9)  - ETA(8))/4 

WD2  = WD2  - PGC(2)  * 40  * (ETA  (6)  4 ETA(24)  4 ETA(25)  4 
ETA(26)  - ETA(8)  - ETA(23)  - ETA(22)  - ETA(2l))/4 
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WD2  = WD2  - PGC(3)  * 4C  *(ETA(21)  + ETA  (22)  + ETA(23)  + 

ETA(8)  - ETA(15)  - ETA(17)  - ETA  (18)  - ETA(19)/4 

WD2  = WD2  - PGC(4)  * 40  * (ETA(8)  + ETA(9)  + ETA(10)  + 

ETA(1)  - ETA  (14)  - ETA(13)  - ETA(14)  - ETA(15))/4 

The  Moment  About  The  Pilot's  Location  is: 

WM6  = PGC{1)  * 40  * (ETA(6)  + ETA(7)  + ETA(8)  - ETA(l)  - ETA(2) 

- ETA(3))  *8/3 

WM6  = WM6  + PGC(2)  * 40  * (ETA(26)  + ETA(25)  + ETA(24)  + ETA(6)  - £TA(21)  - 
ETA(22)  - ETA(23)  - ETA  (8  ))  *50/4 

WM6  = WM6  + PGC(3)  * 40  * (ETA(19)  + ETA(20)  + ETA(21)  - ETA(15)  - 
ETA(16)  - ETA{8))  *28/3 

WM6  = WM6  + PGC(3)  * 40  * (ETA(21)  + ETA(22)  + ETA(23)  + ETA(8)  - 
ETA  (19)  - ETA(18)  - ETA(17)  - ETA(15))  *50/4 

WM6  = WM6  + PGC(4)  * 20  * (ETA(15)  +ETA(16)  +ETA{8)  - ETA(12)  - 
ETA(ll)  - ETA(1  ))  *28/3* 

WM6  = WM6  + PGC(4)  * 40  * (ETA(8)  + ETA(9)  + ETA(10)  + ETA(l  ) - 
ETA  (6)  - ETA<5)  - ETA(4)  - ETA(3))  *40/4 

WM6  = WM6  + PGC{1)  * 40  * (ETA(6)  + ETA(5)  + ETA(4)  + ETA( 3) 

- ETA(8)  - ETA(9)  - ETA(10)  - ETA(l))  *10/4 

WM f'  = WM6  + PGC(2)  * 20  * (ETA{21)  + ETA(27)  + ETA(26)  - ETA(8)  -ETA(7) 
-ETA(6))*  8/3  . 
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Forces  and  Moments  due  to  Cushion  Pressure.  Upward  forces  are  exerted 
on  the  vehicle  at  the  centroid  of  each  cushion  by  the  air  pressure  in  the 
cushion.  The  force  on  each  compartment  is  the  cushion  pressure  (PGC) 
times  the  compartment  area  which  is  800  square  feet.  Pitch  and  roll 
inompntfi  are  also  caused  by  these  forces: 


In  the  esg  direction,  heave: 

4 

FX3CT  = - £ PGC(J)  x 800 

3=1 


(58) 


In  the  pitch  direction,  es-: 

3 


FX4CT  = £ 
3=1 


PGC(J)  x 800  x X.. 


In  the  roll  direction,  es^: 


FX4CT  = £ 
3=1 


-PGC(J)  x 800Yj 


where  Xj,  Yj  is  the  location  of  the  compartment  center  in  plan  relative  to 
the  ship  based  origin. 


Forces  and  Moments  due  to  Spray  and  Skirt  Drag.  The  forces  due  to  spray 
from  the  moving  cushions  and  the  cushion  skirt  drag  are  coupled  in  the 
experimental  work  on  ACVs.  We  have  no  reliable  basis  to  attempt  a sepa- 
ration of  these  two  effects  analytically.  Therefore,  these  effects  have  been 
approximated  by  losses  proportional  to  horizontal  velocity  squared  following 
experimental  data  in  the  Bell  PDSR : 


u 


es, 

j. 

eso 


SDBX  = - CjjXj  lXx| 


(59) 


SDRY  = -CDX2  jX2l 


es. 


3DRR  = -SDRY  (Hp) 


es. 


SDRP  = SDRX  (Hp) 


esi 


SDRYM  = SDRY  (XS1C)  - SDRX  (XS2C). 


The  drag  coefficient  Cq  is  assumed  equal  to  0. 25  based  on  the  Bell  skirt  drag 
god  spray  losses.  Hp  is  the  height  of  the  origin  above  water. 


1 
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ENVIRONMENT  MODELING 

The  components  of  the  vehicle  environment  modeled  are  the  visual 
scene,  wind,  waves,  and  beach.  The  visual  scene  is  created  photographi- 
cally from  a camera  mounted  on  the  gantry.  The  wind  direction  and  velocity 
are  inputs  to  the  program  and  are  combined  with  the  craft  velocity  vector 
to  find  apparent  wind  relative  to  the  craft.  Aerodynamic  forces  due  to  the 
apparent  wind  arc  calculated  by  curve  fits  to  Bell  experimental  results. 

The  wave  surface  is  composed  of  the  seaway,  or  waves  that  were  present 
before  the  craft  arrived,  and  vehicle  generated  waves. 

WIND  AND  AERODYNAMIC  FORCES.  Apparent  Wind.  The  wind  velocity 
in  feet  per  second  (VWIND)  and  the  wind  direction  in  radians  clockwise  from 
North  (AWIND)  are  inputs  to  the  simulation.  From  these  and  the  craft 
velocity  vector  the  wind  relative  to  the  craft  is  calculated  (apparent  wind): 


I O 


Head  wind  component: 

Side  wind  component 
(from  starboard): 

Magnitude  of  apparent 
wind; 


XWIND  * X,  + VWIND  cos  (X_ -AWIND) 
1 b 


YWIND  = X0 -VWIND  sin  (X  -AWIND) 
<5  b 


VA  =A/ XWIND2  + YWIND2 


The  apparent  wind  direction,  positive  clockwise  from  straight  ahead  is 
then  computed: 


For  YWIND  + : 
For  YWIND  - : 


Tan’ 1 (YWIND/ XWIND) 

- 7T+  Tan"1  (YWIND/ XWIND)  . 


f fg 

S fi&Q 


This  angle  in  degrees:  ABETA  3 BETAA  (180/ tr  ) has  a range  of  -180 
to  +180. 

Momentum  Drag.  The  vehicle  experiences  reaction  forces  due  to  the  ingest 
of  air  at  relative  velocity  VA  into  the  cushion  fan  system: 

MDRAG  = P (QFAN1  + QFANg)  VA.  (6- 

This  drag  results  in  forces  and  moments  about  the  pilot  location: 


% £ 
¥ 9^ 

g ?T, 

esl: 

XMDRAG 

3 

-MDRAG  cos  (BETAA) 

(65) 

c 5"* 

®®2: 

YMDRAG 

— 

-MDRAG  sin  (BETAA) 

(66) 

1 ^ 

i %.  a 

•v 

PITCHM 

S 

XMDRAG  (XS3C) 

(67) 
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ea .:  HOLLMD 

4 

es„:  YAWMD 


-YMDRAG  (XS3C)  (88) 

(124XS1C)YMDRAG-XMDRAG(XS2C).(69) 


Windage.  The  side  force  drag  coefficient  i3  curve  fit  from  Bell  experimen- 
tal data  shown  in  figure  3. 

SDRAG  = .01385  |ABETA|-  7.69  x 10"5  ABETA2  (70 

if,  ABETA  < 0,  SDRAG  = - SDRAG. 

The  frontal  drag  coefficient  is  curve  fit  from  Bell  experimental  data  shown 
in  figure  10,  Thi3  curve  is  fit  with  a parabola  for  angles  less  than  4G  , 
with  a cubic  from  40°  to  88°,  and  with  a parabola  above  88  . Pert  star- 
board symmetry  is  assumed: 

if,  iABETA|<40O  (71 

FDRAG  = -2. 22  x 10-4  ABETA24-3. 33xl0"3  |abETa|  4 0. 5 
if  40°<  (ABETA)  < 88° 

FDRAG  = -2,48  x 10~5  [ABETA  | 3 4-  5 x 10'3  ABETA2 
- .324  | ABETA (4-  5.835 
if,  |abETA.|>88° 

FDRAG  * 1.04  x 10~4  ABF,TA2-3.  518x10~2|ABETA|  4-  2.39. 

The  yaw  moment  coefficient  about  the  boat  center  is  based  on  the  Bell 
experimental  data  shown  in  figure  11.  This  curve  is  fit  with  a straight  line 
for  |ABETAj<60°,  with  a cubic  for  jABETA  j between  60°  and  120°,  and 
with  a straight  line  from  120°  to  180°.  Again  port- starboard  symmetry 
is  assumed, 

if,  |ABETA|<  60°:  YDRAG  = 1.67  x 10”3  [ABETA(  (72 

if,  60°  < ABETA<120°:  YDRAG  = 5. 17  x 10‘7|ABETA| 3 - 1.23  x 10"4 

ABETA2  4-  5.82  x 10'3 |aBETa|+  8.27  x 10~2 


ABETA  > 120  : 


YDRAG  * 1.67  x 10 
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The  resulting  windage  forces  and  moments  are: 


es^: 

XBDRAG 

= 

-FDRAG  x F CARE  A x l/2pVa2 

(73) 

**2: 

YBDRAG 

c 

-SDRAG  x SCAREA  x l/2pVa2 

(74) 

^5: 

PITCHA 

= 

XBDRAG  x XS2C 

(75) 

es 

4 

ROLLA 

ss 

-YBDRAG  x XS3C 

(76) 

Tss: 

YAWBD 

YDRAG  x SCAREA  x LCUSH 
x l/2pVa2  + YBDRAG  x XS1C 
-XBDRAG  x XS2C. 

(77) 

Windage  forces  also  occur  on  the  propeller  ducts.  The  angle  of  attack  of  the 
relative  wind  to  the  duct  if  first  calculated  including  yaw  velocity  of  the 
vehicle  and  propeller  slipstream: 


-1  f Va  sin/?aa  +AX  x XD(6)  1 

ADUCTj  = tan  1 Y ^ * 

1 Va  cos  /3aa  + 1. 78  VTROPj 
Duct  sideforce  coefficient: 

CDUCTj  = -7.407  x lO-3  ADUCTj  + 1.333 
if  ADUCTj  <18,  CDUCTj  = 6.667  x 10"2  ADUCTj. 

Effective  velocity  at  duct: 


if.  TPROPj  >0.0 


{va  cos/3aa  + 1.78  VTPROPj  }2  + {va  sin£aaj 


VDUCTj  = 

if  TPROPj  <0.0 

VDUCTj  * y^Va  cos£aa  - 0. 89  \/TPROPj|2  + (va  sin£aaj^ 

Effective  dynamic  head  at  duct: 

VDj  = l/2p  VDUCTj2. 


(78) 


(79) 


(80) 


(81) 
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The  vehicle  forces  and  moments  due  to  the  duct  are  then  calculated: 

YDUCT  = -£2  CDUCTj  x DIA  x CHORD  x VDj  (82) 

ri 

YAWD  = YDUCT  x X1R 

ROLLD  = YDUCT  x X3R  (83) 

The  aerodynamic  yaw  damping  moment  is  estimated  by  Bell: 

YAWDC  = -2. 77  x 10®  x XD(6).  (84) 

SEAWAY  DESCRIPTION.  The  composite  small  amplitude  linear  wave  and 
finite  amplitude  solitary  wave  ocean  surface  model,  as  detailed  in  the  June 
report,  was  subsequently  found  to  be  unworkable  within  a real  time  simula- 
tion environment.  The  essential  difficulty  derived  from  attempts  to  match 
in  both  time  and  space  the  surface  elevations  given  by  the  two  wave  forms. 

A simpler  model,  judged  to  be  adequate  for  engineering  and  simulation 
purposes,  has  been  devised  and  is  detailed  below. 

As  a wave  enters  shallow  water,  roughly  at  a point  where  the  water 
depth  is  one-half  of  the  wave's  length,  the  shoaling  process  begins,  hi 
general,  the  process  reduces  the  wave's  length,  increases  the  height  of  the 
crest,  and  reduces  the  depth  of  the  trough.  Small  amplitude  wave  theory 
can  be  used  directly  to  predict  the  first  effect  within  engineering  accuracy 
and  to  partially  account  for  the  latter  two.  These  affects  are  conveniently 
accounted  for  by  using  the  Stokes  second  order  correction  for  the  wave  profile. 

Offshore  of  the  point  where  shoaling  begins,  the  first  approximation 
to  the  ocean  profile  is 


o H 

V = cos  (K  x - wt).  (85) 

e>  O 

O 

As  the  wave  shoals,  however,  both  the  wave  height  and  the  wave  number  vary 
with  distance,  giving: 


o H (x)  __ 

7)  a — 2 — cog  (K(x)  - wt).  (86) 


The  wave  height  in  shoaling  water  H(x)  can  be  solved  for  in  terms  of  the 
deep  water  wave  height  H , the  local  water  depth  h(x),  and  the  local  wave 
number  K(x).  It  is  (ref.  fppen,  page  65): 


H°(x) 


H 

o 


2 jcob  h^  (K(x)  h(x)  ) 

2 K(x)  h(x)  + sin  h(2K(x)  h(x)  ) 


(87) 
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In  equation  87,  the  local  wave  number  K(x)  is  determined  by  solving  (itera- 
tively) the  transcendental  equation  for  the  dispersion  relation: 

w2  = g K(x)  tan  h (K(x)  h(x)  ).  (8 

The  spatial  variation  as  given  by  the  K(x)  term  in  the  argument  of  the 
cosine  in  equation  86  requires  special  interpretation.  In  deep  water  the 
rate  of  change  of  wave  phase  in  space  is  given  by  K ; thus,  the  phase  is 
given  by  ° 

x 

/K  dx*  * K x (8 

o o 

as  shown  in  equation  85.  For  shoaling_water,  however,  the  analytic  form 
for  this  rate  of  change  is  unknown  and  g(x),  a function  which  must  be  com- 
puted for  each  beach-wave  length  combination,  is  defined  as  follows: 


J5 

Em  ■/ 


K(x»)  dx'. 


As  the  water  depth  decreases,  the  slope  of  K(x)  will  increase  giving  a more 
rapid  spatial  variation  of  the  cosine's  argument  in  equation  86,  equivalent  to 
a shorter  wave  length. 

Equation  85  and  its  inshore  generalization,  equation  86,  gives  a first 
approximation  to  the  wave  height;  one  in  which  the  crests  and  troughs  are 
displaced  an  equal  distance  from  the  mean  surface  elevation.  A second  order 
correction  can  be  made  by  utilizing  results  from  Stokes  second  order  wave 
theory.  This  improvement  is  based  on  the  first  order  solution.  With  crests 
occuring  at  both  the  crests  and  troughs  of  the  first  order  wave,  the  effect 
is  to  increase  the  overall  height  of  the  crest  and  decrease  that  of  the  trough. 
This  second  order  addition  then  may  be  written  • 


rj  = H (x)  cos  [2  (g  (x)  - wt)  J. 
1 


H (x)‘ 
7 r o 


K(x) 


cos  h(K(x)h(x)  )(2+cos  h(2K(x)h(x)  ) 


sin  h (K(x)  h(x)  ) 


To  prevent  this  second  order  contribution  from  erroneously  dominating  the 
first,  H (x)/H°(x)  is  limited  to  a maximum  value  of  1/8.  This  condition 
prevents  additional  crests  and  is  justified  scientifically  by  recalling  that 
this  contribution  is  treated  only  as  a small  correction  term. 
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As  in  the  previous  seaway  description,  the  waves  are  assumed  to 
form  "spilling  breakers"  and  disperse  their  energy  gradually  as  they  approach 
the  shore.  Critical  wave  height  to  water  depth  and  wave  height  to  wave  length 
ratios  are  imposed  as: 


H°(x)  / h(x) 


<0.78 


H (x)  / (2  tt/  K(x)  ) <0.142  . (94 

Again,  there  is  no  wave  remaining  at  the  beach  (by  equation  93)  and; 
therefore,  there  is  no  run  up  to  account  for.  This  departure  from  reality 
should  not  significantly  alter  the  craft  dynamics  or  training  experience. 

Seaway  Calculations.  Calculation  of  the  wind-generated  waves  requires 
executing  an  initialization  subroutine  as  well  as  a sequence  of  calculations 
within  the  time  cycle  loop  in  the  simulation. 

The  main  program  reads  a wave  period,  offshore  wave  height  and  a 
beach  slope,  and  then  calls  the  subroutine  SEAWAY  which  does  the  following: 
Computes 

WLEN  The  offshore  wave  length. 

WNKO  The  offshore  wave  number. 

DZ  The  depth  at  which  wave  shoaling  begins. 

XSHOAL  The  distance  offshore  at  which  shoaling  beging 

(measured  value  of  the  X01  coordinate). 


L * 


The  remaining  portion  of  the  program  determines  values  of  the  function  K(x) 
as  described  by  equation  90  in  Seaway  Description.  This  is  done  using  a 
crude  rectangular  integrating  routine  with  step  sizes  in  the  S direction  equal 
to  the  deep  water  wave  length.  (As  written  only  100  values  are  permitted; 
this  could  be  increased  if  necessary).  During  execution  of  SEAWAY,  the 
subroutine  WVNBR  (for  wavenumber)  is  frequently  called.  This  routine 
solves  the  transcendental  equation, 

w2  = gK(x)  tan  h lK(x)  h(x)  j 

relating  the  finite  depth  water  wave  frequency  to  its  wave  number  K(x)  and 
the  water  depth  h(x).  This  is  done  using  a Newton-Raphson  root  seeking 
technique  for  the  equation  (given  w and  h) 

F(k)  = w2  - g K tan  h (Kh).  (95 

Once  SEAWAY  has  been  executed,  no  further  calculations  need  be  made 
prior  to  the  beginning  of  the  time  loop. 

Within  the  subroutine  OCEAN,  the  wind  wave  heights  are  calculated 
at  35  pts  and  saved  in  the  array  ETA  (I) ; I = 1,  35. 
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Each  of  these  points  correspond  to  a point  of  water  surface  directly 
below  a known  point  on  the  craft.  • Since  the  beach  has  been  assumed  to  run 
east-west  (X02  direction),  the  land  to  occupy  X-l  >0,  and  the  waves  to 
approach  the  beach  perpendicularly,  the  wave  shape  depends  only  upon  the 
X01  coordinate  and  time.  When  given  a point  In  space,  the  first  step  is 
to  make  certain  that  it  is  over  water,  not  on  the  land. 

If  (HPXOl(N)  ) 10,  50,  50 

1 If  on  land,  HPX01  >0  and  ETA  (N)  is  set  to  zero. 

If  at  sea,  a further  check  is  made  to  see  if  the  point  is  in  the  wave  shoaling 
region  or  further  offshore. 

If  offshore,  the  wave  height  and  length  are  unaffected  by  the  beach  and  ARGU; 
the  argument  of  the  cosing  giving  the  shape  of  the  wave  is  simply  calculated. 

If  in  the  shoaling  region,  the  wave  height  is  adjusted  according  to  the  Seaway 
Description  and  the  spatial  contribution  to  the  argument  determined  by  linear 
interpolation  within  tne  table  generated  by  subroutine  SEAWAY. 

Finally,  the  Stokes  Second  Order  corrections  are  made  on  either  the  offshore 
or  shoaling  'waves  again  according  to  SEAWAY  description.  This  process 
is  repeated  for  each  of  the  35  field  points. 

VEHICLE-GENERATED  WAVES.  The  prediction  of  the  surface  wave  pattern 
generated  by  the  motion  of  the  craft  is  a complex  problem.  The  complexity 
lies  in  the  fact  that  the  craft's  motion  is  not  known  beforehand,  and  the  tra- 
jectory of  motion  in  the  horizontal  plane  could  be  very  arbitrary.  Difficult 
as  the  problem  is,  nevertheless,  it  is  necessary  to  come  up  with  a wave 

. height  production,  because  vehicle-generated  waves  not  only  provide  the 

basic  ingredient  for  calculating  wave  drag  and  moment,  but  also  carry  the 
effects  of  changing  the  cushion  volume  as  well  as  the  skirt  clearance,  fac- 
tors which  crucially  affect  the  response  of  the  craft. 

The  mathematical  model  for  the  prediction  of  vehicle-generated 
waves  is  based  cm  linearized  water-wave  theory.  The  formulation  and 
derivation  of  the  solution  of  this  hydrodynamics  problem  are  detailed  in 
Appendix  A,  The  main  result  is  that  to  obtain  the  wave  height  at  a given 
point  wider  the  craft,  a convolution  integral  in  time  must  be  evaluated. 

Input  to  the  evaluation  of  the  integrand  and;  thus,  the  integral  consists  of 
the  recenv.  history  of  the  trajectory  and  orientation  of  the  craft  as  well  as 
the  coordinates  of  tne  point  where  the  wave  height  is  desired.  The  integrand, 
or  the  kernel  f inaction  of  the  integral,  is  a aery  complicated  mathematical 
expression,  evaluation  of  which  during  real  time  simulation  is  certainly 
beyond  the  present  state  of  the  art  of  computer  technology.  To  make  real 
time  evaluation  of  the  convolution  integral  possible,  tabulation  of  the  kernel 
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function,  which  is  a function  of  three  variables,  is  necessary.  With  the 
table  of  kernel  function  loaded  in  memory,  the  computer  can  be  program- 
m ed  to  interpolate  values  from  the  table  and  evaluate  the  integral  on  a 
real  time  basis. 


Therefore,  to  implement  the  proposed  mathematical  model  for  the 
prediction  of  vehicle -generated  waves,  we  can  divide  the  required  program- 
ming efforts  into  two  stages,  hi  Stage  one,  values  of  the  kernel  function 
are  generated.  This  is  an  expensive  and  time-consuming  process.  The 
mathematical  technique  and  associated  programming  procedures  are  tho- 
roughly documented  in  Appendix  B.  These  calculations,  however,  have  to 
be  done  only  once  for  a given  water  depth.  A medium  size  table  has  been 
generated,  for  straight-ahead  motion  as  well  as  a limited  amount  of  sway 
and  yaw.  This  table  can  be  extended  with  little  conceptual  difficulty  to 
accommodate  the  case  of  a general  maneuver. 


i 


hi  Stage  two,  which  occurs  during  the  real  time  simulation  process, 
the  convolution  integral  which  depends  on  the  trajectory-history  of  the  craft 
as  input  will  be  performed.  The  numerical  procedure  as  well  as  the  program- 
ming efforts  involved  In  this  stage  are  described  in  considerable  detail  in 
the  main  body  of  this  report.  Examples  and  test  results  are  given  to  illus- 
trate the  application  of  the  kernel  table  generated  in  Stage  one.  The  presen- 
tation that  is  to  follow  consists  of  two  parts:  The  first  part  describes  the 
principles  and  mathematical  details  of  the  numerical  procedure  used  in 
evaluating  the  convolution  integral  with  specific  reference  to  the  origin  of 
the  equations.  The  second  section  consists  simply  of  the  FORTRAN  IV 
version  of  the  procedure  described  in  the  first  section,  put  into  effect. 

Results  for  several  test  runs  of  constant  and  variable  speed  conditions 
are  given.  Finally,  to  interface  the  vehicle-generated  wave  package  with 
the  rest  of  the  motion  simulation  program,  we  must  provide  a routine  that 
transfers  and  updates  the' array  which  stores  the  history  of  the  craft's 
trajectory  and  orientation.  This  is  described  in  Section  HL 


f ' 
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NUMERICAL  PROCEDURE  FOR  REAL-TIME  COMPUTATION  OF  VEHICLE  - 
GENERATED  WAVES.  The  nondimens ional  vehicle-generated  wave-height, 
't'*  f(x’,  y')/(p  / pg)  can  be  obtained  by  performing  the  following  convolution 
integral  with  respect  to  the  nondimens  ional  time,  *x  = t / Yajg  : 

a ° 

f(x’,y’)  = f dr  KS(P,  (r),  R (r),  r*).  (96) 

x y 

o 

where  (x’,y‘)  are  the  coordinates  in  the  horizontal  plane  at  which  the  wave 
height  is  desired,  and  Ks  is  the  kernel  function  which  will  be  tabulated 
and  stored  in  a common  block  described  later.  To  evaluate  (96)  we  must 
also  know  the  history  of  motion  of  the  craft  from  t *=  0 (i.  e. , the  present) 
to  r = -T  , say  some  NT  units  of  time  interval,  AT,  back  in  the  past.  This 
is  reflected  througirthe  functional  dependence  of  Ry( T ) and  Rx(?^  ) on  1*, 
as  given  by  equation  (7.4)  of  appendix  A: 

Rx(^  ) = Ix*  - X(  )]  cos6(^  ) + [y‘  - Y(  )]  sin  6 (r  ), 

Ry(^r  ) = Iy‘  - Y(^)J  cos  ${  ^)  - [x’  - X(^)]  si ).  (97) 

Here,  X(f ),  Y(7>),  and6(T)  are  the  history  of  the  craft’s  trajectory  and 
orientation,  measured  in  the  craft’s  present  coordinate  system,  and  must 
be  known  before  we  can  calculate  the  ^-integral.  All  length  quantities 
shown  have  been  nondimens ionalized  by  the  craft’s  half-length,  or  . 

0>f  00  00 

Let's  assume  we  know  X{  "f  ),  Y(  r ),  6 { r ) in  a discrete  manner 
as  follows: 


fl 
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Here,  we  use  the  two-dimensional  array  STATE  (n,  m)  to  store  the  history 
of  the  trajectory.  Note  that  if  we  integrate  NT  time  intervals  of  AT  each 
in  the  past,  there  are  actually  NT+3  ordinates  because  we  use  half-intervals 
in  the  recent  past.  Note  also,  by  definition,  STATE  (o,  n)  for  m*l,2,3. 

Assume jthat  we  are  given  x*,  y:  and  the  STATE'S  are  defined,  then 
the  integrand  K&(R  (t*),  R (t*),  t)  at  , n=0,  ...  , NT+2,  can  be 
generated  by  calling  a function  subprogram  K SIN TP  (RX,  RY,  N): 

REAL  KSINTP,  F(50) 

NN  = NT  + 2 
DO  10  N*l,  NN 

^Calculate  RX,  RY  according  to  Equation  (97) 

I using  STATE  (N,M),  M=l,2, 3 and  x',  y'.  (99) 

AX  = ABS  (RX),  AY  = ABS  (RY) 

10  F(N)  = KSINTP  (AX,  AY,  N) 

The  integrand  is  now  stored  in  the  array  F;  the  wave  height,  £ , can  now 
be  obtained  simply  by  applying  Simpson's  rule  with  half-interval  adjustments 
as  below: 

?(x\y')  = -~{  | F(0)  + 2*  (F(l)  + F(3)  ) + F(2)  +|f(4)  + F(NN) 

+ 4*[F(5)  + F(7)  + F(9)  + . . . + F(NN-l)  ] (100) 

+ 2*{F(5)  + F(8)  + F(10)  + . . . + F(NN-2)  jj 
where  F(0)  = 0,  and  NN  must  be  even. 

We  now  proceed  to  describe  the  interpolation  details  of  the  function 
subprogram  KSINTP,  assuming  tha*  a table  of  the  kernel  function  is  available. 

KSINTP  (AX,  AY,  N)  Routine.  The  Kernel  function  KS(x,y,  ir)  is  a function  of 
three  variables.  One  may  envision  this  as  a three-dimensional  array,  with 
subscripts  i,  j,  n,  corresponding  to  grid  points  indices  of  the  variables  x, 
y,  and  ^respectively.  The  table  is  prepared  in  such  a way  that  for  given 
values  of  Rx,  r , and  r,  only  interpolations  in  the  x and  y directions  are 
necessary.  In  other  words,  the  t*  grid  of  the  Ks  table  is  identical  to  the 
discrete  time  values  defined  in  (98). 

A schematic  diagram  showing  the  storage  pattern  of  the  values  of 
the  kernel  function  is  shown  on  the  following  page: 
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In  order  to  completely  define  the  table,  the  following  information 
should  be  stored  in  a common  block  called/KRNL,/  : 


V.. 

ijn 


value  of  the  kernel  at  the  i,  j grid  pojnt  on  the  table 

corresponding  to  the  time  value  T - r . 

n 


(NX) 


(NY) 


x and  y grid  point  values  of  the  -table. 


Number  of  x-grid  points  in  the  -table. 
Number  of  y-grid  points  in  the  -table. 


For  example: 

COMMON/  KRNL/  NX(5  0),  NY  (50),  X(32,  50),  V(32,  15,  50) 


Each  call  of  K SIN  TP  (RX,  RY,  N)  will  therefore  interpolate  the 
value  of  the  Kernel  function  from  the  x and  y grid  points  in  the  -table. 
The  value  is  returned  through  the  function's  name.  The  following  inter- 
polation algorithm  is  based  on  a 3-point  Lagrange  formula. 
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Step  1;  Search  and  find  the  indices  p and  q such  that 

X ^ RX  < X , andY  <RY<YJ, 
p,n  p+l,n  q,n  q+l,n» 


If  p = NX  (N)  - 1, 

If  q = NY  <N)  - 1, 


set  p = p-1* 
set  q = q-1. 


Step  2:  Interpolate  the  value  of  kernel  for  3 y-grid  values  at  the  x-value 
we  want. 


t— r RX-X.  rx  RX  - X. 

G e f]  l ( HL)}*v  + [*l  { — )]*  V 

j VXin  P’q4,n  V.H-V 

r-r  RX-X. 

+ 1 ii^i.  (x—  Tor*  1 * for <10i> 

p p ’ 

Step  3;  Interpolate  in  the  y-direction  to  obtain  the  Kernel  value: 


.1.  q+2  RY  - Y. 

KSINTP  ■ 2_, 1 IT  < Y * * Gk! 

k=o  q+k-1'  jn 

3#(q+k) 


(102) 


Note  that  considerable  sin?  'lifications  in  these  formula  will  occur  when 
the  grid  spacings  are  coi.  mt. 

Subroutines  for  the  Computation  of  Vehicle-Generated  Waves,  The  "vehicle - 
generated  wave  package n which  contains  four  subroutines  is  basically  a 
FORTRAN  IV  version  of  the  procedure  described  in  the  first  section  above. 
The  name  and  purpose  of  eac'  of  the  components  inside  the  package  are 
given  below: 

Routine  LDTABL.  This  routine  should  be  called  once  at  the  beginning  of  the 
calling  program  which  uses  the  package.  Its  purpose  is  simply  to  load  the 
common  Dlock/KR'NIi/  by  reading  the  deck  of  cards  which  contains  the 
values  of  the  kernel  function  Ks. 

Routine  VGWAVE.  This  subroutine  performs  the  convolution  integral  in 
time  to  obtain  the  wave  height.  Before  entry,  the  STATE  array  must  be 
defined  f >r  at  least  as  many  time  steps  as  the  kernel  table  is  prepared. 
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In  this  routine,  the  procedure  given  by  Equations  (9S)  and  (100)  of  the  section 
above  is  performed  by  calling  the  interpolation  routine  KSXNTP  to  generate 
the  integrand  and  the  routine  STMPSN  to  carry  out  the  integration. 

Routine  KSINTP.  A function  subprogram  which  interpolates  a value  of 

the  kernel  function  from  the  kernel  table  stored  in  the  common  block  /KRNL/. 

Routine  SIMPSN.  General  purpose  Simpson's  rule  integrator. 

The  above  routines  are  given  together  with  a sample  main  program 
which  calls  VGWAVE  to  obtain  the  wave  Height.  In  this  particular  applica- 
tion, the  STATE  array  was  generated  by  assuming  that  the  craft  is  going  at 
constant  speed.  Typical  results  are  shown  in  figure  12  and  figure  13  for  a 
number  of  speeds.  The  Froude  number  corresponding  to  the  case  where  the 
depth  Froude  number  equals  1.0  is  0 • 707,  above  which  waves  diminish  in 
height  substantially.  It  has  been  observed  that,  in  doing  these  computations, 
a nondimensional  time  step,  A T , of  0 » 25  is  too  large  for  speed  with  Froude 
number  greater  than  5,  although  such  step  size  appears  to  be  more  than 
sufficient  for  low  spee«  cases.  The  kernel  table  used  in  these  calculations 
was  generated  up  to  T =6.0  (see  Equation  (96))  back  in  the  past.  This  is 
insufficient  for  speed  lower  than  F =0*5,  while  for  high  speed,  a value  of 
T = 2.0  is  sufficient.  For  an  accuracy  of  ar,roximately  5%,  a useful  cri- 
terion for  the  truncation  time  TQ  of  the  convolution  integral  (Equation  (96)  ) 
is  to  keep  integrating  until  [r|(T0)  + R|(T0)  ] a 25.0.  Physically,  this 
implies  that  effects  caused  by  the  craft  more  than  2. 5 craft-lengths  ago 
are  neglected.  Such  capability  can  be  easily  incorporated  in  the  final  simu- 
lation program. 

One  interesting  test  computation  performed  as  a check  of  the  results 
was  the  case  of  zero  forward  speed.  The  results  for  integrating  up  to  T =6. 0 
are  shown  in  figure  14.  This  resembles  and  approaches  the  shape  of  the  static 
wave  profile,  which  is  entirely  consistent  with  the  fact  that  if  pressure  im- 
pulses are  applied  continuously  over  the  same  spot  on  the  free  surface,  the 
transients  will  eventually  die  out  and  one  is  then  left  with  just  the  static  profile. 


KIDSHIP  C CUT  M-PK3 
rraft  at  Constant  Speed 


(B,  PIST  fM  MIDSHIP  fFR.  NO.  a 
Figure  14. : Craft  at  Zero  Forward  Speed 
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sample  program  for  using  vehicle-generated  waves  package. 

»»«**»•««••  PVY  8/4/?4c 


PACKAGE  CONTENTS 


ROUTINES  »VGWAVE*»  *LDTABL*,  *KTNSTP*»  *SIMPSN*. 
AND  DECK  OF  CARDS  CONTAINING  KERNEL  VALUED. 


THE  FOLLOWING  COMMON  - STATEMENTS  SHOULD  BE  INCLUDED  IN  CALLING 
PROGRAM. 

COMMON/RNPARM/ALPP*STfDFLT.NTSfEPR*  HT 
COMMON /STATUS/ AT  *BT .TIMF .STATF (50,3) 

COMMON  /KRNL/  NX (30) » NY(30),  X(40,30),  Y(?0»30)»  V(40,20,30) 
1 NTAU*  Y (30) 

DIMENSION  Ft.DPT(50.2>  , HEIGHT  (50)  » FTaNON<50> 

OEFINE  A NUMBFO  OF  CGNSTANTS. 

THESE  CONSTANTS  ARE  NOT  ACTUALLY  USED,  BUT  SERVE  TO  DOCUMENT  T«E 
GENERAL  ENVIRONMENT  IN  WHICH  THE  KERNEL  TABLE  WAS  PREPARED. 
ALPP=DECAY  FACTOR  IN  PRESSURE  DISTRIBUTION, 

ST  = Bt'AM/LENGTH  RATIO 

HT  = DEPTH  OF  WATER  / HALF-LENGTH  RATIO. 

ERR  “ ERROR  BOUND  IN  GENERATING  TABLE. 

ALPP=  .314159 
ST  = .5 
ERR  = .01 
HT  a 1. 


OEEINE  PHYSICAL  CONSTANTS  OE  WA 
PRES?  = NOMINAL  PRESSURE  OF  CUS 
RHO  1.98 
GRAVTY  = 32.17 
PRESZ  = 109.375 
CONVRT  = ORES7/(RH0»GRAVTY) 


CONSTANTS  OE  WATER  AND  CRAFT. 
PRESSURE  OF  CUSHION,  IN  LBS/ (SO.  FT) 


CALL  *LDTA0L*  ONCE  TO  LOAD  CARD  DECK  CONTAINING  TABLE  . 

SEE  ROUTINE  *LDT ABL°  FOR  DFTAIL  DEFINITION  OF  VARIABLES  IN  THE 
COMMON  /KRNL/. 

THE  TIME  GRID  WILL  BE  BE  PRINTED  OUT  AUTOMATICALLY. 

IPRINT  SHOULD  BE  SET  EQUAL  TO  1 IF  ONE  WANTS  TO  PRINT  OUT  THF 
VALUE  OF  THE  ENTIRE  TARLF,  OTHERWISE,  SET  IPRINT=0. 

IPRINT  = 0 

CALL  LDTABL  (IPRINT) 

OELT  = T(NTAU)  -T  (NTAIJ-I ) 


GENERATE  SOME  TYPICAL  FIELD  POINTS  DATA. 

NFP  = NO.  OF  FIELD  DOTNTS. 

ARRAY  *FLOPT*  CONTAINS  THE  X AND  Y COORDINATE^  OF  THE  POINT  WHERE 

THE  WAVE  HEIGHT  IS  DESIRED. 

NEP  = 50 

DX  = .1 

BY  = .5 

00  100  J=1 *2 

DO  100  K=l,  2S 

KK  = U-l>*25  ♦ K 

FL0PT(KK,1)  = 1.  -(K-1)*0X 


2001 


FL0PT(KK,1)  = 1.  -(K-1)*0X 
FLOPT (KK .2)  = (J-1)*DY 
CONTINUE 

WRITE (6,2001 ) ( (FLOPT(K,J) , J=l,2),  K,  K=I,  NFP) 

FOPMAT (//*  FIELD  POINTS,  X AND  Y »/  (2X,  2Fln.4,l5)) 

DO  400  M=l,  10 


GENERATE  SOME  TRAJECTORY  DATA  AND  STORE  IN  ARRAY  «STATE« 

IN  THE  FOLLOWING  TEST-DATA  GENERATION,  CONSTANT  SPEED  IS  ASSUMFD. 
NTS  = NO.  OF  TIME  STEPS  FOP  WHICH  STATE  DATA  ARE  AVAILABLE  , AND 
MUST  BE  .SE.  TO  NTAU 
READ  (5,3001)  FRNO 
IF  < FRNO  .LT.  .01  ) STOP 
3001  FORMAT  ( F10.4) 
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) 


WRITE  (6,  300?)  FRNO 
3002  FORMAT  { *1  FROUDE  NO.  = ••  F5.3 

SQ2  - SORT (2.0) 

NTS  " NT  AU 

NTS 

STATE (N,l ) = -SQ2*FRN0*T(N) 

STATE (N.2)  = 0. 

200  STATE (N.3)  a 0. 

WRITE (6*2092)  { (STATE ( I * J) * J=1 *3) « 1=1, NTS) 

2002  FORMAT  { *0  STATES  W.R.T.  PRESFNT  COORDINATE-SYSTEM*,/ 
1 * XOR  YOR  THPTA  * /,  (1X.3F10.5  )) 


«*• 


WITH  FLDPT  ANO  STATE  DATA  DEFINED.  WE  ARE  READY  TO  CALL  *VGWAVE* 
TO  OBTAIN  WAVE-HEIGHTS. 


SAMPLE  LOOP  FOR  USING  ROUTINE  *V GWAVE*.  WHICH  RETURNS  A NON- 
OIMENSIONALIZED  WAVE-HEIGHT  VALUE. 

MULTIPLYING  BY  THE  VARIAELF  *C0NVRT*  WIIL  CONVERT  RESULT  TO 
PHYSICAL  UNITS. 


»««•# 

IPRINT  = 0 
DO  300  1=1, 


NFP 


IPRINT.  IF  SET  EQUAL  TO  1.  WILL  CAUSE  AIL  INTERPOLATED  KERNEL-VALUES 
TO  BE  PRINTED  OUT  WHEN  *VGWAVE»  IS  CALLED.  IF  THAT  IS  NOT  DESIRED. 
SET  IPRINT=0. 

ETANON(I)  =-VGWAVE(  FL0PT(I,1),  FLDPT C I *2).  IPRINT) 

HEIGHT (I)  = CONVRT*ETANON(!) 


300  WRITE  (6,3004) 
«««** 


I,  FLDPT (1,1),  FLDPT(I,2),  HEIGHT(I),  ETANON(I) 


3004  FORMAT  ( *0  FIELD  POINT  NO.*, 13,*  (X,Y)=<»,  2F8.4, * ) • ,3X , 


1 • WAVE  HEIGHT  =•*  F*.4,?X,  * LENGTH  UNITS*, 

2 *.  (N0N-DTMENSI0NALI7ED=*,  F10. 6,  «)*  ) 


STATEMENT  BELOW  SERVES  MERELY  TO  PUNCH  RESULTS  ON  CAROS. 

WRITE  (7,  3005)  FRNO,  (I,  FLDPT (1,1),  FLDPT (I, 2).  ETANON(I) , 
1 1=1,  NFP) 

3005  FORMAT  ( *FP.  NO.  =*,  F«.4,  / (15,  5X*  2F10.4,  F10.6)) 

400  CONTINUE 


STOP 

END 


m b 
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C 

C 

C 

C 

C 

C 

C 


FUNCTION  VGWAVFt  XFP*  YFP,  IPRTNT) 

COMMON/RNP ARM/ALPP * ST . DEL  T » NTS  * FRR * HT 
COMMON/STATUS/AT. BT,TIMF,ST4IF<50, 3) 

COMMON  /KRNL/  NX<30).  NY(30)»  XC40.30)*  Y<20*30)»  V<40,20,30). 

1 NTAU*  T (30) 

REAL  KSINTP 
DIMENSION  F (50  ■ 

THIS  SUBROUTINF  PERFORMS  THE  CONVOLUTION  INTEGRAL  IN  TIME  TO 
OBTAIN  THE  WAVE-HFIGHT  AT  THE  FIELD  POINT  WITH  COORDINATES  < XFP*  YFP) 
BEFORE  ENTRY*  THE  ARRAY  *STATE*  DESCRIBING  THE  TRAJECTORY*  MUST  RF 
DEFINED  FOR  AT  LEAST  *NTAU*  TIME  STEPS. 

THE1 COMMON  BLOCK  /KRNL/  TS  ALSO  ASSUMED  TO  HAVE  BEEN  LOADED  ALRtADV. 
IF  ( IPRINT  .EO.  1 ) WRTTF  (6*  1002) 


• J«,  RX*  toX(TAU)**  AX*  *RY(TAU)  »*6X»  »TAU*»  4X* 


NTAU 


1002  FORMAT ( //3X  * 

1 • KERNEL*  ) 

F ( 1 ) = 0. 

DO  30  N=l* 

COMPUTE  THE  QUANTITIES.  PX(TAU)*  AND  RY(TAU). 

DT  = STATE (N. 3) 

CD  = COS(OT) 

SD  = SIN(OT) 

DX  = XFP  - STATE (N,l) 

DY  = YFP  - STATE <N, 2) 

AT  = CD*DX  ♦ SD*DY 
BT  = CD*DY  - SD*DX 
RX  = A8S ( AT) 

RY  = ABS(BT) 

GENERATE  KERNEL  VALUES  BY  CALLING  KSINTP. 

JJ  = N* 1 

F ( JJ)  = KSTNTPfRX,  RY,  N) 

TF  ( IPRImT  .NE.  1 ) GO  TO  30 
WRITE  (6.  1003)  N*  AT,  BT,  T (N) , F(JJ) 

1003  FORMAT  ( IX,  14,  2X,2F11.3,  2X,  F8.4  , 2X*  F10.6  ) 
30  CONTINUE 


CALL  INTEGRAION  ROUTINE 

CALL  SIMPSN  (F,  CELT , NTAU*1* 

VGWAVE  = ANS 

RETURN 

END 


ANS) 
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REAL  FUNCTION  KSINTP<RX*  RY « N) 

INTERPOLATION  ROUTINE  FOR  KERNEL  VALUE.  _ 

COMMON  /KRNL/  NX(30)»  NY (30)*  X <40*30)*  Y<20*30)*  V(40*20*30)* 
1 NTAU*  THO) 

INTEGER  P«  Q 

STEP1..  SEARCH  FOR  THE  INDICES  P AND  0. 


JI  = NX  CN) 
0 10  1=3 


GO  TO  25 


DO  10  1=3.  II 

P=I 

IF  ( X(PtN)  .GE.  RX)  GO  TO  15 
10  CONTINUE' 

GO  TO  50 
15  Pap-2 

JJ  * NY (N) 

DO  20  J=3.  JJ 

0= J 

IF  ( Y (Q«N)  .GF.  RY)  GO  TO  25 
20  CONTINUE 
GO  TO  50 
25  0=0-2 

STEP2*  INTERPOLATE  VALUES  OF  3 Y-GRID  VALUES  AT  THE  VALUE  OF  X GIVFN. 

RP  a RX-X(P.N) 

RPla  RX-X(E»*1«N) 

RP2  a RX-X(P*2.N) 

PP1  a X(P«N)-X<P*1.N) 

PP2  = . X (P»N) -X (P*2tN) 

P12  a X(P*1 *N)-X  (P*2*N) 

Fl=  RP1*RP?/{PP1*PP2) 

F2=  -RP*PP2/(PPi*P12) 

F3  = RP*RP1/IPP2*P12) 

GO  a F1*V(P*Q  *N)  ♦ F2*V(P*1*0  *N)  ♦ F3*V(P*2*Q  *N) 

G1  = F1*V(P.QM  *N)  ♦ F2*V(P*i*G*l.N)  *_?3*V CP*2*Q*1 *N) 

G2  = F1*V  CP *0*2* N) ♦ F2*V (P«-l  .0*2*N)  ♦ F3*V(P*2.Q*2*N> 

STEP3  * CROSS  INTERPOLATE  IN  THE  Y-DIRECTION. 

RO  = RY-Y(O.N) 

R01  a RY-Y(0*1.N) 

ROE  a RY  - Y (Q*2*N) 

001  a Y(0*N)  - Y(0*1.N) 

002  a YIQ»N)  -Y(Q*2*N) 

012  a Y(Q*1.N)-Y<Q*2*N>  


KSlNTPa  R01*RQ2/ (001*002) *G0 
L RO*R01/ (002*012) *G2 


- R0*RQ2/ (001*012) *G1  ♦ 


RETURN 
50  KSINTP=0. 
RETURN 
END 
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SUBROUTINE  LDTABL (IPRINT) 

: ROUTINE  FOR  LOADING  KERNEL  TABLES. 

COMMON  /KRNL/  NX(30).  NY(30),  X(40*30>*  Y<20.30)»  V<40,20,30). 
1 NTAU.  T(30) 

T (N)  = TIME  VAl  UE  CORRESPONDS  TO  THE  N-TH  TANLE. 

NTAU  = NO.  OF  TIME  TABLES. 

X = X-GRIO  POINT  VALUES.  Y = Y-G»IO  POINT  VALUES. 

NX  CN)  = NO.  OF  X-GRIO  POINTS  IN  TABLE  N. 

NY <N)  = NO.  OF  Y/GRID  POINTS  IN  TABLE  N. 

V a KERNEL  VALUE. 

1001  FORMAT  l IS) 

1002  FORMAT  (25X,  F10.5) 

1003  FORMAT  ( 32X,  15  /(8F10.5)) 

1004  FORMAT  (8X«  9F8.4) 

1012  FORMAT  <//•  KERNEL  TABLE  FOR  TIME  = '.  F10.4  ) 

1013  FORMAT  ( » RX  GRID.  NO.  OF  GRID  POINTS  =»,  15.  /,  32X.  IS/ 

1 (8F10.5  ))  ^ 

1023  FORMAT  { » RY  GRID.  NO.  OF  GRID  POINTS  = *.  15.  /.  32X,  I*;/ 

1 (BF10.5  )) 

% 

READ  (5.1001)  NTAU 
DO  10  N=l.  NTAU 
READ  (5.  1002)  T (N) 

11  READ  (5.  1003)  II.  (X(I.N).  1=1.  II) 

12  NX (N)  = II 

READ  (5.  1003)  JJ.  (V(J,N).  J=l.  JJ> 

13  NY (N)  s jj 

IF  ( IPRINT  .NF.  1 ) GO  TO  15 
WRITE (6.  1012)  TCN) 

WRITE (6.  1013)  II.  (X ( I »N) * 1=1.  II) 

WRITE  (6.  1023)  JJ.  (Y(J.N).  J=l.  JJ) 

15  DO  20  J=l.  JJ 

READ  (5.  1004)  ( V(I.J.N).  1=1,  II) 

IF  ( IPRINT  .NF.  1 ) GO  TO  20 
WRITE (6*  1004)  ( V(I.J.M).  1=1.  II) 

20  CONTINUE 
10  CONTINUE 

WRITE  (6.  1005)  NTAU,  ( T (N) * N=l,  NTAU) 

1005  FORMAT  ( 2X*  15.  * TIME  GRID  VALUES  AS  DEFINED  BELOW..*, 

1 /,  C8F10.4  )) 

RETURN 

END 
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v-Cfi  SUBROUTINE  SIMPSN  <Y*  H*  N*  R) 

REAL  Y ( 1 ) • SE*SO*  R 

M=(N-3)/2 

so=o. 

v-  S£=  Y (2) 

fF  (N  .EO.  3)  60  TO  15 
00  10  1=1*  M 

SO  = SO  ♦ Y {2*T*1 ) 
w 10  SE  = SE  ♦ Y(2*I*2> 

15  R=  H* ( YC1)  *Y(N>  ♦ 4,*SE  ♦ 2.*S0  )/3. 
RETURN 
END 
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The  Interface  Between  Vehicle  Kinematics  and  the  Computation  of  Vehicle- 
Generated  Waves.  Determination,  of  the  vehicle-generated  wave  system 
requires  a knowledge  of  the  past  path  arid  orientation  of  the  vehicle.  This 
section  outlines  two  mathematical  manipulations  required  to  provide  the 
vehicle-generated  wave  package  the  necessary  information  regarding  the 
vehicled  trajectory:  (1)  conversion  to  the  coordinate  system  used  in  the 
vehicle  generated  wave  subroutines  and  updating  the  records  with  most  re- 
cent motions;  and  (2)  providing  this  information  in  dimensionless  form 
and  at  time  increments  as  needed  by  the  vehicle -generated  wave  subroutines. 

The  first  operation  is  accomplished  by  the  subroutine  UPDATE  whose 
argument  list  contains  newly  determined  incremental  changes  of  the  vehicle’s 
position  and  orientation  in  the  horizontal  plane,  as  measured  in  ihe  coordinate 
system  fixed  at  the  pilot's  location.  This  routine  updates  the  dimensional 
records  of  the  vehicle  motion  which  is  maintained  in  the  coordinate  system 
used  in  the  analysis  of  vehicle  generated  waves.  The  pilot  based  system 
has  XI  forward,  X2  starboard  and  X6  (yaw  angle)  measured  positive  from 
the  XI  axis  toward  the  X2  axis;  the  system  used  in  the  determination  of 
vehicle-generated  waves  has  X forward,  Y port  and  6 yaw  angle  measured 
positive  from  the  X axis  toward  the  Y axis  (figure  15a).  Because  the  coor- 
dinate system  in  which  the  motion  records  are  maintained  is  continuously 
changing  (fixed  to  moving  craft)  the  entire  record  needs  to  be  updated  at  each 
time  step.  This  can  be  seen  clearly  in  (figure  15h). 

During  a time  interval  between  time  t and  time  (t-Mt)  the  integrated 
equations  of  motion  will  give  an  incremental  change  of  XI,  X2,  and  X6, 
say  DX1,  DX2  and  DX6.  From  these  subroutines  UPDATE  first  determines 
corresponding  changes  DX,  DY,  DD  according  to: 

DX  = DX1  + 30.0  - 34.9857  cos  (DX6  + 0.54042026) 

DY  = -DX2  - 18.0  + 34.9857  sin  (DX6  + 0.54042026} 

DD  = -DX6 


) 


v_ 


The  equations  reflect  geometrical  relations  best  explained  by  refer- 
ence to  figure  15a.  The  next  operation  refigures  the  motion  history  in  the 
new  coordinate  system  located  DX  ahead,  DY  to  port  and  rotated  DD  to  port 
from  the  old  system.  This  record  is  kept  in  the  array  DSTATE  (100.3) 
where  the  columns,  left  to  right,  contain  the  x position,  y position  and  yaw 
angle  of  the  vehicle  at  a number  of  time  intervals  into  the  past  corresponding 
to  the  row  number.  (Note  that  row  "0"  isn't  needed  since  the  three  values 
there  are  by  definition  zero. ) During  each  time  step,  the  coordinates  of 
a previous  point  on  the  trajectory  makes  one  step  down  the  column.  The 
oldest  entry  is  always  lost.  Using  expressions  for  the  shift  and  rotation 
of  a coordinate  system's  origin: 
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The  second  operation  is  performed  by  the  subroutine  INTERP.  This 
routine  uses  array  ESTATE  to  produce  array  STATE,  which  differs  in  two 
ways.  First,  the  entries  are  made  in  "natural"  length  units  "a”  (vehicle 
Half  length)  and  in  radians.  Second,  the  time  increment  is  changed  to 
correspond  to  the  dimensionless  time  used  in  the  vehicle-generated  w ave 
package  and  the  time  increment  between  all  entries  is  not  the  same.  Be- 
tween entries  in  rows  one  through  four  the  time  increment  is  1/8  of  a 
"natural"  time  unit  ■/a/g  (a  = vehicle  half  length;  g «=  gravitational  acceler- 
ation). Between  rows  four  and  five  (and  thereon)  the  tune  increment  is  of 
a time  unit. 
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SECTION  IV 

LOGIC  FLOW  FOR  ACV  MATHEMATICAL  MODEL 


Initialize  for  Stopped  Craft 


Set  up  seaway  characteristics 


Get  wavej  *}.,  i=  1,  35  for 
seaway  from  OCEAN 


Relative  Wind 

Eqn  60 

, 61,  62 

Get  Prop  Torques  and  Powers 
Eqn  13,  14 


Get  £ an  Torques  and  Powers 
Eqn  29,  30 
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Calculate  cushion  pressures 
and  flows 

Eqn  48  to  52 


Calculate  vehicle  generated 
wave  heights  & sum  with  seaway 

Eqn  54 


Momentum  Drag 
j Eqn  64  to  69 


Wind  Forces 

j 

Eqn 

70  to  77  | 

— — » 

Duct  Forces 

Eqn  78  to  83 
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APPENDIX  A 

SOLUTION  TO  THE  HYDRODYNAMIC  PROBLEM  OF  A PRESSURE 
DISTRIBUTION  MOVING  WITH  ARBITRARY  SPEED 
ALONG  AN  ARBITRARY  PATH 

The  following  summarizes  all  the  important  mathematical  details  on  the 
formulation  and  solution  of  the  problem  of  a pressure  distribution  moving  in 
any  arbitrary  manner  in  the  horizontal  plane.  The  craft  is  allowed  to  have 
any  speed  and  any  heading  angle.  The  assumptions  are  that: 

a.  The  fluid  is  inviscid 

b.  The  flow  is  irrotational  and  potential;  and, 

. c.  The  waves  generated  by  the  motion  of  craft  is  small  enough  so  that 

the  linearized  free-surf ace  condition  can  be  used. 

1.  Coordinate  Systems  and  Other  Definitions 

x,  y,  z are  coordinates  based  on  inertial  frame  of  reference 

x‘,  y’,  z‘  are  coordinates  referred  to  a moving  system  which  is  mounted 
on  the  craft.  Location  and  orientation  of  this  system  is  given 
by  X(t),  Y(t),  m 

z axes  point  upwards 


O 


Trajectory  is  given  parametrically  as  function  of  time,  t: 
t t 

X(t)  = / U(t)  d t = / [C  (r)  cos  8(r)  - C (r)  sin  8]  dr, 
0 0 x 7 

t t 

Y(t)  = f V(t)  dr=  f [C  (r)  sin  8(r)  + C (t)  cos  8]  dT, 
0 0 x 7 

Cx  = U cos  8 + V sin  8 , 

C = - U sin  8 + V cos8  . 

y 


(103) 

(104) 

(105) 

(106) 
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Relations  between  the  two  systems  are: 

£=  £'cos  8(7)  - V sin  8<r)  + X(r). 

^ = |*  sin  8(r)  + cos  $(r)  + Y(t). 

and 

£'  = ft-  X(t)]  cos  8 + h-  Y(t)]  sin  8, 

V1  = -[  |-X(t)]  sinS  + [■>?-  Y(t)  ] cos  8 - 

2.  Formulation  in  the  Inertial  Coordinate  System 

The  pressure  distribution  in  the  moving  coordinates  is  assumed  to  be 
given  as  p'(£  ij  ’,  r)  in  general,  in  the  inertial  system,  pft,  rj , t).  The  region 
in  which  of  p (£,17  , t)  is  nonzero,  changes  with  time.  The  fluid  is  assumed  to 
be  of  a constant  depth  h,  and  density  P.  Define  the  velocity  potential 
= $(x,  y,  z,  t),  then  continuity  equation  gives  y2$=  0. 

Boundary  Conditions 

On  Free  Surface  F:  (z  = 0) 

$tt(x,  y,  0,  t)  +g$z(x,  y,0,  t>-pt(x,  y,  t)/P.  (109) 

On  Bottom,  fi:  (z  = -h) 

$2<x,  y,  -h,  t)  = 0.  (110) 

Here,  g is  the  acceleration  of  gravity. 

These  conditions  occur  very  frequently  and  derivations  are  available  from 
any  text  book  on  free-surface  hydrodynamics.  See  for  instance,  J.  J.  Stoker 
(1957,  p.  191). 

Initial  Conditions 


(107) 


(108) 


$t(x,  y,  0,  0)  = -g#x,  y,  0)  - p/p  (x,  y,  0)  = 0,  (111) 

$z(x,  y,  0,  0)  = £t(x,  y,  0)  = 0.  (112) 

I.  E.  No  initial' disturbances  except  a static  deformation  of  the  free  surface. 
Here  £ is  the  wave-height. 

3,  Solution  in  terms  of  the  Time-Dependent  Green  Function 
As  usual,  one  starts  with  Green’s  Theorem  for 

4ttA(x,  y,  -z,  t)  = jJ[G(x,  y,  z;  £,  *7,  £ ; t,  t)0  -<p  G,  I dS.  (113) 

1 *"fubuzr  v 1 
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where 


SL  - ( JjL,  } * n;  G is  a Green  Function  with  a 1/r 

singularity,  and  Sr  is  a large  vertical  cylinder  enclosing  all  the  disturbances. 

It  can  be  shown  that  the  integral  over  £„  vanishes  in  the  limit  as  R,  the  radius 
of  the  cylinder,  approaches  infinity,  ana  after  a number  of  operations*  on  (113) . 
we  may  write: 

4»r#x,y,  z,t)  = J d7Gt(x,y,z,£,^,0,T,t)ptU,J?,T).  (114) 

-00  ” 

In  arriving  at  (114),  we  have  made  use  of  the  inital  and  boundary  conditions 
of£  and  have  assumed  that  G satisfies  the  following  properties: 

G = 1/r  + H(x, y, z;  £,??,£,  t ,r), 

V2G  = 0,  for  (x,  y,  z)  *(£,»?,{);  V2H  = 0, 

G*'iG‘‘U=0' 

G^(x,y,z;  V,  -h;  t,r)  = 0, 

G(x,y, z;  q,£,  t,t)  = Gt(x,y,z;  £,*?,;,t,t)  = 0.  (115) 

Basically,  G satisfies  the  free-surface,  the  bottom  boundary  conditions,  and 
two  initial  conditions  when  the  source  is  brought  into  existence  at  time  equal 
to  r.  The  Green  function  was  first  derived  by  Lunde  (1951),  [ by  Haskind 
(1946)  for  the  infinite  depth  case]  and  is  available  from  Wehausen  & 

Laitone  (1960,  p.  494). 


In  the  present  notation,  the  Gr  sen  Function  takes  the  fo  llowing  form: 
G(x,  y.  z;  4.  4;  t,  T0)  = 1/r  * l/r2  - 2 J ">*  ' 

cosh  k(z  + h)  J0  (kR)  dk  + 2 f dk^ik  — 'MgL,  f dr 


sin  [7(k)  (t-T)]  cosh  k(h+£)  Jg(kR) 


(116) 


~see  for  instance,  J.  J.  Stoker,  Water  Waves,  1.  9,  p.  187, 
Inter  science  Publisher,  1957. 
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where 

r2  = (x  )2  + (y-V)2  -?•  (z-£r, 
r|  = (x  - £)2  + (y-vT  + (z  + 4)w, 

r2  = (x  -o2  + (y  -y)2, 
y(k)  = ygk  tan  h kh 

and  J„  is  the  Bessel  function  of  order  zero.  For  our  application,  we  are 
concerned  only  with  £ = 0,  and  G therefore  simplifies  to 

f oo 

G(x,  y,  z;  y,  0;  t,  tq\  = 2 J dk  C0^^kf^h)  [1  - cos  y{k)  (t  -rQ)] 


JQ  (kR)  dk. 


(117; 


Differentiating  this  and  substituting  it  into  equation  (114),  we  get  the 
following  equation  for  <£(x,  y,  z,  t): 


5>(x,'y,z,t)  J dr  jj  d£cP)p  (^,t)  J dk 

&rr  n ct-r-s  * n 


q -S(T) 

cosh  k(z4h)  sin'/  \t-~Q 
cosh  kh  y 


JQ(kr) 


(118) 


4.  Solution  in  the  Moving  Coordinate  System 

Since  p{£,  1),  r)  can  be  expressed  as  pM^1,  y')  in  the  moving  coordinates, 
in  which  the  integrals  o\  er  (|,  h)  do  not  have  limits  depending  on  time,  it  is 
more  convenient  to  write  (118)  in  terms  of  x1,  y',  z’.  Using  the  coordinate 
transformation  developed  earlier,  and  noting  that  in  the  moving  coordinate 
system: 


d?  _ dP'x  dP’ST+lP’^V 
Sr  dr  ' £7}  ^ * 


(119) 


it  can  be  shown  that  after  some  fairly  straight  forward  calculations: 

o-oo  y ch  kz 


dusinm-r)  chk(_z+h)  ft  d , 
y ch  kz  ^ 


- i p*(  77’,t)  *[w(U<T)  - (^' sinS +>?'cos8)5)  + u(V(t) 

i ar 


+ (f’cosS-^1  sin8)S))p  ciw  (X(x,,y,,t)  - X (£*,  r)] 

+ iu  lYfc’.y'.t)  - Y(Z',y',T)]}  i = 


(120) 
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where 

X(x’,y*,t)  = X'  cosS(t)  - Y*  sinS(t)  + X(t) 
Y(x*,y',t)  = X*  sinS(t)  + Y*  cosS(t)  + Y(t). 


\121) 


In  arriving  at  (120),  we  have  used  an  integral  representation  of  Jn(kR) 
to  convert  the  k-integral  to  the  w-u  double  integral.  Here  U 

2 A 2 .2 
w + u = k . 

Equation  (120)  is  the  solution  to  the  problem  of  an  arbitrarily  moving 
and  spinning  craft.  This  is  the  analog  of  Equation  (2-21)  of  the  report  by 
Doctors  and  Sharma  (1970),  now  extended  for  3 degres  of  freedom  instead 
of  just  rectilinear  motion. 

Note,  if  the  pressure  in  the  moving  system  is  nortime  dependent,  the 
first  term  drops  cut.  Note  also  if  S =0,  S * 0,  this  reduces  to  a solution 
corresponding  to  purely  translational  motion  with  a drift  angle  but  no 
spinning.  If  we  had  formulated  the  problem  in  the  moving  coordinates. 
Equation  (120)  will  fce  the  solution.  A good  check  will  be  to  start  with  (120) 
and  see  if  the  free-surface  condition  in  the  rotating  system  will  be  satisfied. 
The  free-surface  condition  is  the  version  linearized  with  respect  to  dis- 
turbance caused  by  p;  no  linearization  with  respect  to  the  smallness  of  angular 
rotation  is  necessary. 

5.  General  Expression  for  Wave-Height 

To  obtain  the  wave-height,  one  simply  uses  equation  (111)  now,  [RHS^Oj 
which  relates  the  wave-height  to  the  time -derivative  of  the  potential,  but 
must  bear  in  mind  the  meaning  of  the  time  differentiation  in  a moving  frame. 
Upon  performing  an  integration  by  parts  of  the  4E  term,  with  respect  tor,  we 
obtain  the  following  simplified  expression:  cfr 

«*’•»*••>  'erg  — ^w  du||  d£' d7?1  P,(£',1?,,t>  cosXt  Exp{  i\v[X(x'y't)  -£*J 

+ iutYfc'.y'.t)  -V']\ — R dw  du JJ  df  d*j’  t dr 

4J7pg  -00  S'  0 

p(C,  rj\r)y.  sinr(t-T)  * Exp{iw  [X(x',y\t)  - X<£’, ^*,r)l 
+ iu  [Y(x«,y*,t)  - Y(£,,7?,,r)]f.  (122) 

with  X(x,y,t)  and  Y(x',y’,t)  defined  as  before. 
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TMs  expression  has  the  advantage  that  one  only  needs  to, keep  track  of  X(t), 
Y(t),  andg(t),  instead  of,  in  addition,  Cx(t),  C (t),  and  S(t).  Equation  (122) 
i3  the  wave-height  seen  in  the  craft's  frame,  \rave  generated  as  a result  of  a 
given  history  of  motion  prescribed  by  X(t),  Y(t),  and  S(t). 

To  evaluate  some  of  these  integrals  analytically,  we  must  now  know  the 
functional  Jorm  of  pHVtVW- 

6.  Functional  Form  of  the  Pressure  Distribution 


Equation  (122)  i«  a quintuple  integral.  The  (C',1?')  integration  corresponds 
essentially  to  taking  the  two-dimensional  Fourier  transform  in  space  of  the 
pressure-distribution  function  p (£,*‘>7,*t).  The  (w,u)  integrals  correspond  to 
the  inversion  of  the  transforms,  and  the  time  integral  represents  an  inte- 
gration of  the  effects  created  by  the  craft's  trajectory  and  orientation  in  the 
past.  It  should  be  noted  that  for  the  .simulation  problem  being  considered  at 
hand, ' it  is  impossible  to  calculate  the  time  integral  ahead  of  time.  This  is 
so  because  the  trajectory  and  orientation  of  the  craft  depends  on  the  actions 
taken  by  the  pilot,  actions  which  are  unknown  before  the  simulation. 

i 

Because  of  the  immense  generality  of  the  problem  and  complexity  of  the 
integrals,  it  will  be  to  our  advantage  to  be  able  to  perform  the  (£',■>?')  integrals 
analytically.  We  suggest  here  adopting  the  hyperbolic -tangent  distribution 
used  by  Doctors  and  Sharma  (1970).  The  functional  form  of  p^'.'n'.r)  is 
given  by 

p(C.  V'.t)  * ?G  [tanh«(  I'  + 1)  - tanh  ot(  |*  - 1)] 

x [tanh  P ( | ) - tanh  £(  J*  - | )]  (123) 


where  a,  p are  diffusion  parameters  in  the  longitudinal  and  the  transverse 
directions  respectively;  the  larger  they  are,  the  closer  the  shape  is  to  that 
of  a step  distribution,  a,  b in  (123)  are  the  half-lengths  and  half  beams  of  the 
craft  respectively.  P(t)  is  a time  modulation  factor  which  can  be  and  will 
later  be  set  equal  to  1. 0 if  the  pressure  is  assumed  to  be  nontime  dependent 
in  the  craft's  frame,  and  p^  is  the  nominal  pressure.  For  this  distribution, 
pQ  is  given  by 

P()=A/4ab  (124) 

where  A is  the  displacement  of  the  craft. 

With  this  assumption  of  the  functional  form  of  p(£',i? ’t)  the  (£',  Y) 
integrals  can  be  evaluated  analytically,  as  done  earlier  by  Doctors  and 
Sharma  (1970,  pp.  37-40).  The  result  is 
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I \ p(  >7*  )e*  cos&*  u sinS)^  (-w  sinS  + u cos  8)1?  ] dV  d^' 

-TO-OO 


ir  P0P<T)  sin  aw'  • sin  bu* 


<*P 


sin  h(jif~)  sin  h (^) 


(125) 


U 


where  the  notations 

w*  = w cos  8 + u sing, 

u*  = -w  sin8  + u cos  8 , (126) 


have  been  used,  bearing  mind  w and  u are  variables  of  integration  in  the  next 
step. 


Next,  substituting  (125)  back  into  (122)  and  changing  variables  of  integration 
from  (w,  u)  to  (w*,  u*)  as  defined  in  (126)  we  get 

£ = *P<°>  KC  *1 

- J*  dTP(T)  KS  [R  (t-T),  R (t-T),  t-r] , (127) 

o x y 

where  the  Kernel  Junctions  are  defined  as  follows: 


K 

K' 


s(Rx-  Ry  »* 


IS 

>00 


dw'  du'  e1 


+ iu'R  / cos  y t\  sintew'jsintu'b) 
^ resin  y tj 


«0sh  (^wl)Sh(^u*) 


with  (128) 

y = Vgk  tanh  kb  , k2  = w'2  + u'2  • 

C c 

It  can  be  seen  that  K and  K are  functions  of  three  variables  for  given  values 
of  a and  /3.  Also,  the  quantities  R (t-T),  R (t-r)  are  merely  symbolic  since 
they  do  not  depend  explicitly  on  thexcombinafion  (t-T).  Physic  ally,  they  have 
the  meaning  of  being  the  x-  and  y-  coordinates  of  the  point  (x*,y')  in  the  craft's 
system  at  time  equals  t,  expressed  in  the  coordinate  system  of  the  craft  at 
time. equals,  to  t.  For  a given  set  of(x',y')  at  time  t,  R (t_T)  and  Rv(t~r) 
may  then  be  calculated  as  follows:  x 

Ex  = [X(x',y«,t)  - X(t)J  cosS(t)  + lY(x’,y*,t)  - Y(r)j  sin  8(r), 

R = - [X-X(r)j  sin  S(t)  -!•  lY-Y(r)]  cos  B(t), 

v 

where  (x,y)  are  given  by  Equations  (121). 
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It  is  worthwhile  to  point  out  at  this  point  that  a good  choice  of  the  values 
for  « and/3  are  a=0=5.  0.  As  discussed  in  Doctor's  and  Sharma's  report,  this 
choice  though  appeared  arbitrary,  provide  two  advantages.  First,  it  effec- 
tively eliminates  the  unrealistic  hollows  and  humps  of  the  wave  resistance  curve 
at  low  speeds.  Second,  the  computational  efforts  involved  in  evaluating 
the  wave-number  integrals  can  be  reduced  considerably  because,  as  can  be 
seen  from  Equation  (128),  the  decay  factor  is  dependent  on  and  ipn  for 
the  (w-u)  integrals. 

7.  Nondimensional  Expression  for  the  Free-Surface  Elevation 


Equation  (127),  developed  earlier  for  the  free-surface  elevation,  consists 
of  two  terms.  The  term  associated  with  the  function  can  be  thought  of 
as  the  wave  system  generated  as  a result  of  the  initial  building  up  of  the 
cushion  pressure.  This  system  of  waves  can  be  neglected  if  the  forward 
motion  of  the  craft  does  not  occur  "immediately”  after  the  build  up.  If  the 
first  term  is  neglected,  what  is  left  then  is  a convolution  integral  in  time 
involving  the  Kernel  function  Ks.  This  integral  can  be  evaluated  ii  the 
quantities  X(t),  Y(t),  and  5(t)  are  given. 

In  actual  computations,  it  is  much  more  convenient  to  use  a nondimen- 
sional form  of  Equation  (127)  ard  measure  time  i with  the  present  being 
equal  to  zero.  Therefore  nondimensionalizing  all  physical  lengths  by  the 
craft’s  half-length  a,  and  time  by^a/g,  for  example, 

*x  = x'/a 

^ t 

W — W * P 

f = (t-fl/vOF  !/, 

*y  - (5 tank  1?  (h/a)»J  ' , etc.,  (129) 

(127)  becomes: 

X&yX-  0)  = s(x',y',t)/(p0/4pg) 

=:]  drPff)  KS(Rv,  R . r)  (130) 


where 


K8^,  Ry,  t) 


*»  . <w  ,b. 
sin  w smu  f— j 


7 sin  rt- 


a0sh(fa)sh^). 


(131) 


with 


aRx(r)  = [X’-X(T)j  cos  S(r)  + [Y«-Y(T)j  sin  S(r), 
aR  (?)  = - [X'-X(r)]  sin  8(r)  + [Y'-Y(r)]  cos  8(T). 


(132) 
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In  writing  down  (132),  we  have  assumed  that  X(?),  Y(r)r  8(r*},  are 
given  in  the  craft's  present  coordinate  system.  In  other  words,  X<r-0)  - 
Y(¥=0)  = $(t  =0)  = 0,  X(^)  and  Y(r*)  are  negative  if  the  craft  is  moving 
forward  and  swaying  to  port;  and 8 (*r)  is  positive  if  the  craft  is  turning  to 
starboard.  In  so  doing,  we  can  dispense  with  the  inertial  coordinate  sys- 
tem used  earlier  in  formulating  the  problem  and  it  is  then  only  necessary 
to  know  the  history  of  the  trajectory  and  orientation,  measured  relative 
to  the  present  position,  rather  than  the  absolute  value  measured  in  an 
inertial  frame. 
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APPENDIX  B 


ti 


THE  KERNEL  FUNCTION,  K‘  (R  ,R  ,t) 

x y 

In  this  appendix,  mathematical  and  numerical  techniques  for  eval- 
uating the  kernel  function  K (R  , R ,t)  are  described.  Consider  the  double 
integral  defined  by  eqn  (131)  of’app^enciix  A. 


oo 

KS(R  »Rv»t)  = j'j' dwdu  g^wI*x+uRy^  • Y sin  yt  • 


sin  w • sin  — u 
a 

a/)sh  (¥F  > sb  <2j 1 


I'.*. 

K (R 

X 

■f- 

1; 

where 

Here  for  simplicity  of  notation,  we  omit  all  the  titles  above  the  dimensionless 
variables  with  the  understanding  that  all  quantities  under  consideration  in  this 
appendix  are  dimensionless. 

THE  KS  INTEGRAL  IN  POLAR  COORDINATES.  Assuming  0 = 5,  i.  e. , the 
diffusion  parameter  is  the  same  in  both  longitudinal  and  transverse  directions 
and  writing  the  (w,  u)  integral  in  polar  coordinates  (k,  6 ),  we  have 

JP°  cos(kR  cosfi)sin(kcos9)cos(kR  sin0)sin(k/ASin6?) 

,t)-  4 /kdkrsinyt  / d 9 £ 2 W 


a 


b/a 

7r/2a 


a sin  h (k'cosQ)  sinh  (k'sint?) 

beam  to  length  ratio 

a decay  factor,  characteristic  of 
the  distribution 


(133) 

(134) 


k‘ 

a'k  . 

we  define: 

R13 

= [(R  +1)2+(R  - fj-  )2J 
A y 

1/2 

» 

j = 1,2 

R2j 

■ HR.  -i>2+<R  Wj 

" J 

1/2 

* 

3 = 1.2 

+ 

t. 

-1'  Ry  ” ^ 

= tan  (-*■—■)  , 

* = [0., 

13 

+ 

R - 1 

X 

(135) 


R..  and  'f'..  are,  basically,  polar  coordinates  of  (R  , R ) referred  to  the 
U 13  s x y 

four  corners  of  the  craft.  K can  now  be  written  as  follows: 
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;■ 


i ; 

I 


KS(R  ,R  ,t> 


oo 


- o 


sinyt 


{-)1+jQ(k,R  , f..)J 


where 


Q(k,  R,  ) 


1_  cos  {kRcos  (6 -ft)]  - cos  [kRcos  (g-t-ft)  j 

2 ^ sin  h (k5  cos  Q)  sin  h (k*  sing) 


(136) 


(137) 


Our  immediate  task,  evidently,  is  to  be  able  to  compute  the  theta-integral 
defined  by  (137). 


THE  THETA-INTEGRAL  Q(k,  R , ft ).  It  can  be  seen  from  (137)  that  the 
Theta-Integral  is  a highly  oscillatory  integral,  particularly  for  large  values 
of  kR.  The  magnitude  of  the  integrand,  however,  is  exponentially  small  for 
large  values  of  k’.  To  avoid  handling,  this  oscillatory  integral  numerically, 
it  was  found  that  the  following  Bessel  series  expansion  yields  2xcellent 
accuracy  with  minimal  amount  of  numerical  computations: 

oo 


Q(k.  R,ft  ) *2 


Ji 


<-)“  C (k1)  J_  (kB)  • sin  2nft 

n*"l|  b.  • * 


(138) 


where  J2n  is  the  Bessel  function  of  the  first  kind,  of  order  2n,  and 


cn(k') 


sin  2n6dg 

sin  h (k'cosg ) sin  h (k’sinfi)  * 


k'  t 0 . 


(139) 


This  coefficient,  in  general,  has  to  be  evaluated  numerically,  but  is  a function 
of  two  parameters  only,  n,  k*.  The  Bessel  functions  can  be  computed  by 
standard  mathematical  library-routines. 

Because  n is  odd  in  (138),  one  sees  immediately  that  Q(k,R,ft  ) is  an 
even  function  of  about  ft  -tt/  4;  5tt/4;  9ir/4,  and  llr/ 4,  and  odd  about  ft  = 0, 
r/2;  t;  3tt/4.  For  example,  we  expect  the  behavior  of  Q as  a function  f 
as  follows: 


The  behavior  of  Q in  [0,  tt/ 4]  defines  its  entire  behavior  completely.  The 
number  of  terms  that  one  needs  to  take  in  (138)  will  depend  on  the  magnitude 
of  kR.  To  estimate  the  number  of  terms  required,  we  note  that  for  a given 
value  of  z,  a large -n  approximation  for  Jn(z)  is: 


* 
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j„<«>  = -4-  t-f-  >n 

n , rr~  -sn 
V 2-rrn 


Hence,  n can  be  so  chosen  that 


1 GZ 

2 log  (2ir n)  + n log  { ) < log  <■ 


where  f is  a measure  of  the  truncation  error  for  terminating  the  series 
with  a finite  number  of  terms,  e.g. , * = 10~6.  It  is  also  evident  from  (140) 
that  n in  (138)  should  be  greater  than  (ekR  / 4). 

Analysis  of  the  coefficient  C (k1)  indicates  that  Q'k,  R,  9 ) behaves 
approximately  as  follows: 

Q(k-R'*)=  k^SThk'  n4i.3.5...Cn(k,,J2n(kK»Sin(2”*>'  (1 

«!• 

where  C^(k’)  is  "almost"  a constant  function  of  k’.  Note  that  Q decays  like 

l/k’ek*.  The  advantage  of  (138)  lies  in  that  the  oscillatory  behavior  of  Q 
is  kR  is  completely  contained  in  the  Bessel  functions. 

* It  may  be  shown  that  C (k‘)  has  the  following  asymptotic  form  for 
large  value  of  n:  1 

©o 


S?«.Cn<k’)  = FskhF  f*3"11  (F->  ■ ***' 


ir  . nir  V e~s^JL  . 

T,smTs=i T_ A 


which  indicates  that  Cn(k’>  approaches  a constant.  Such 

behavior  is  to  be  expected  because  this  integral  is  reminiscent  of  the  more 

well-known  integral 


f 4 sinnS. 
J sin  2 Q 


d8=  — , for  all  n = odd. 


It  was  found  that  this  asymptotic  formula  is  exceedingly  good  for  n > 10.  In 
practice,  for  low  values  of  n,  Simpson’s  rule  was  applied  to  evaluate  C (k*) 
numerically  over  a range  of  values  of  k’,  (k’=  [.001,  125.]  ).  For  each" 
value  of  n,  the  function  Cp(k’)  is  fitted  with  a spline  curve,  and  the  spline 
coefficients  are  stored  on  cards.  This  idea  results  in  a substantial  saving 
of  computation  time.  For  large  values  of  n,  (n  > 10)  the  formula  (143)  is 
used.  Only  the  first  term  in  the  approximation  was  necessary.  Having 
determined  the  number  of  terms  required  in  the  series,  and  having  gene- 
rated the  CR(k’),  n = 1,3,5,...  either  by  interpolation,  or  by  using  (143), 
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we  reduce  the  calculation  of  Q(k,  R,  'f')  to  merely  a summing  process 
as  defined  by  (138). 

The  routine  package,  Q24,  described  above  puts  into  practice  the 
technique  described  herein.  It  also  takes  advantage  of  the  fact  in  the  double 
sum  of  (136): 


^ Y <-)i+jQ(k,R.,,^.)«X  C (k*)  . iy  y J^(kR..)sin(2n^.)).  (145) 

iSl  Vl  3 3 n-1,3,5. .. n 


2 2 


Thus  redundant  calculations  of  C (k*)  can  be  avoided. 

n 


NUMERICAL  COMPUTATION  OF  THE  K-INTEGRAL.  To  obtain  the  values 
of  K®,  it  is  necessary  to  integrate  the  Q-function  with  respect  to  k.  In 
practice,  the  infinite  upper  limit  has  to  be  replaced  by  a finite  value.  Since 
Q(k,R,t{0,  according  to  equation  (142),  behaves  approximately  as  [k’e^1]”^ 
for  large  k',  we  may  estimate  the  truncation  error  e , incurred  as  a conse- 
quence of  stopping  the  integration  at  kQ,  by  the  following  formula: 


£52  •C*y|sinn  |k^iP=<l)2(¥->1/2/,^re'k'dk'- 

O 


2 2 1/2  - k 

- (— ) k ' e 2a 

” o 


(146) 


For  a value  « « 10"^,  k amounts  to  31. 8,  which  yields  a maximum  value  of 
kQR  = 750,  if  R is  less  ^han  or  equal  to  23.5.  The  value  kQR  is  important 
because  it  determines  the  number  of  terms  required  in  the  Bessel  series 
defined  by  equation  (138). 

The  integrand  of  the  k-integral  defies  simple  analytical  description. 
Oscillation  arises  from  both  the  sin  at  term  as  well  as  from  the  function 
Q(kR).  It  was  observed  that  the  period  of  oscillation  of  the  Q-function  is 
approximately  2ir , with  respect  to  the  variable  kR.  This  is  merely  an 
approximation  because  this  function  is  hardly  sinusoidal.  In  order  to  mini- 
mize the  number  of  evaluations  of  the  integrand,  which  requires  a consider- 
able amount  of  computational  efforts,  it  is  necessary  to  have  some  know- 
ledge of  the  approximate  separation  between,  say,  successive  peaks.  The 
following  model  was  found  useful  in  subdividing  the  interval  k=lO.  ,k  J into 
a number  of  subintervals.-  Let  K*  satisfy  the  relation:  u 

k*  / tanh  k*  (h/a)  = 1.5 

then  for  k<  k*,  an  approximate  "period"  of  oscillation  pr  is  given  by 
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p =2  tt  ! [R*  + VE  tj, 
r 

* . 2 2 ,l/2  , 

where  R = IRX  + Ry  J , h is  the  water-depth/  half-length  ratio,  and 

t is  the  nondimensional  time  parameter  in  the  kernel.  For  k > k 
Pr  = 2 tt  / Ir"  +t  jifk  J 


(148) 


(149) 


In  each  of  the  subintervals,  determined  either  by  (148)  or  (149),  the 
Gauss-Legendre  quadrature  formula  is  applied.  This  was  found  to  be  the 
most  efficient  means  of  performing  the  k-integral  numerically,  without 
having  to  have  to  sacrifice  accuracy.  For  an  absolute  accuracy  of  10~4 
quadrature  formula  of  order  10  for  k < k*  and  of  order  7 for  k<  k*  were 
found  to  be  acceptable  for  this  purpose. 


TABULATION  OF  THE  KERNEL  FUNCTION,  K (R^,  Ry,  t).  In  the  dis- 
cussions above,  we  have  outlined  the  numerical  procedures  for  evaluating  KS. 
The  computations  are  so  involved  that  real  time  evaluation  of  Ks  is  beyond  the 
question.  The  approach  we  suggest  is  to  calculate  ana  tabulate  this  function 
before  the  time  of  simulation.  This  table,  which  is  three-riimensic..al,  can 
be  stored  in  memory  and  then  an  interpolation  procedure,  as  detailed  in 
Section  IV,  can  be  used  to  obtain  the  proper  value  of  K . For  convenience, 
we  recall  the  definition  of  K : 


kH-V>  ‘ -5T 


£x>  2 2 

Jdk  k/sin/t  [ ^ 
o 1=1 


(-)1+jQ(kR..,a k,  .)]  , 
x3  x3 


(150) 


where 


oo 


Q(kR,k',+  ) 


* - o c (k»)  J_  (kR)  sin  2n  f , 
n if  v#  v«  m u 


(151) 


y(k)  = Vk  tanhk(h/a)  , (152) 

and  Rj.,^ij,  are  given  in  terms  of  Rx,  Ry  by  equation  (135)  above.  The 
physical  meaning  of  Ks  is  as  follows:  If  an  impulse  pressure-distribution 
of  the  shape  given  below  is  applied  on  the  free  surface  at  time  equal  to  zero 
in  the  form  of  a delta  function,  KS(RX,  Ry,t)  describes  the  response  at  the 
time  t,  of  the  free  surface  located  at  (R  , R ).  The  free-surface  elevation, 

C,  is  given  by  -(p0/pg)Ks.  Therefore,  K£  is  proportional  to  the  wave  height 
and  contains  all  the  response  characteristics  of  a fluid  with  a constant  but 
finite  depth.  From  this  interpretation,  one  sees  immediately  that  the  region 
of  nonzero  disturbance  will  grow  in  time.  According  to  linearized  water- 
wave  theory,  we  know  that  the  wave  .front,  the  wave  with  the  longest  wave- 
length, will  propagate  the  fastest  and  that  the  maximum  speed  is  at  most 
Vgh,  where  h is  the  water  depth. 


u 
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^ l/2 

In  nondimensional  form,  this  speed  is:  V = (h/a)  ' . Thus,  for  a given  value 
of  t,  one  may  estimate  the  size  of  the  region  in  which  Ks  is  nonzero  as  follows: 


R = Vt 
o 


(w 


If  the  pressure  distribution  had  been  a step  function  (xc,  yc}=(1.0,  p ).  For 
the  distribution  that  we  have  used,  a good  choice  is  {4, , 2, 5). 


The  program  TABKS,  a listing  of  which  is  attached,  was  developed 
for  the  purpose  of  generating  the  Kernel  function  Ks.  Referring  to  (150), 
we  see  that  the  important  physical  parameters  are: 

h * water  depth  to  half-length  ratio 

<*  = diffusion  parameter  of  the  distribution 

H-  = beam/length  ratio  of  craft 

In  addition  to  these  quantities,  we  must  specify  the  array  of  values  of 
RX,R^  and  t,  for  which  the  function  Ks  is  to  be  calculated.  The  computer 
program  TABKS  ready  in  the  array  for  t,  TMS,  the  array  for  Rx,  RS,  and 
the  array  for  Ry  RY.  k , defined  by  equation  (146)  is  then  determined  by  ' 
calling  the  routine  FJNrfko.  Next,  for  each  combination  of  the  values  of 
Rx,  Ry  and  t,  routine  INTDK2  is  called  to  evaluate  the  k- integral  according 
to  the  stens  described  above.  The  results,  i.  e. , the  values  of  Ks,  are  then 
punched  cut  on  cards  (File  #7).  Two  options  are  available  for  using  INJDK2, 
If  the  parameter  1ACY  is  set  equal  to  1,  an  accuracy  of  approximately  -10?e 
will  be  used  as  criterion  for  calculating  the  k- integrals.  This  option-  is  use- 
ful when  some  rough  ideas  about  the  behavior  of  K5  is  desired.  If  IACY  is 
set  equal  to  2,  accurate  calculations  of  more  than  1%  accura  cy  will  be  per- 
formed, and  € in  (146)  will  be  set  equal  to  10-4  for  the  determination  of  kQ. 
Evidently,  in  order  to  obtain  the  integrand  for  the  k-integral,  the  theta- 
integral,  Q(k,  R,»JO  must  be  evaluated.  This  is  achieved  by  calling  the 
routine  Q24  (see  the  fourth  to  last  statement  in  routine  WNING),  which  is  a 
major  component  of  the  Q24  package  detailed  in  the  last  section  of  this  ap- 
pendix. On  an  IBM-370  model  65  machine,  typical  computation  time  for 
each  value  of  Ks  ranges  from  0. 5 to  5. 0 seconds  depending  on  the  magni- 
tude of  the  values  of  R and  t (see  equation  (148)  ). 


hi  figures  B^l  tc  B-13,  we  have  plotted  graphically  the  behavior  of 
Ka  vs.  R , for  different  values  of  R . Each  graph  corresponds  to  a dif- 
ferent value  of  the  nondimensional  time.  These  values  were  generated  using: 


= l.C; 


W 10; 


0.5 
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which  is  applicable  to  the  ACV  being  considered.  The  hump  speed  corresponds 
to  this  particular  geometrical  configuration  as  approximately  15.0  knots.  The 
table  was  generated  for  time-values  of  t=0. 125,  0.25,  0.375,  0.5,  0.75,  1.0, 

(0. 25),  6 0.  In  the  notation  of  Section  IV,  AT  is  0. 25.  The  range  of  Rx 
values  was  chosen  to  include  the  entire  region  of  nonzero  disturbance.  How- 
ever, the  range  of  R is  only  [1. , 1.5],  Hence,  this  table  is  not  sufficient 
for  the  application  or  the  general  case  of  a radical  maneuver,  which  requires 
R to  have  the  same  range  as  F^.  However,  a small  amount  of  s\  -ay  and 
yaw  is  permissible. 

The  following  points  will  be  useful  for  future  preparation  of  extension 
of  the  kernel  table: 

a.  A grid  spacing  of  Rx,  or  R of  0. 1 is  quite  sufficient  if  linear 
instead  of  parabolic  interpolation  is  used  in  the  procedure  described  in 
Section  IV.  More  sparse  grid-spacing  can  be  used  if  higher  order  inter- 
polation function  is  used.  This  becomes  a tradeoff  problem  between  the 
storage  requirement  and  the  available  computational  time  for  performing 
the  convolution  integral  during  real  time  simulation. 

b.  The  value  of  the  kernel  function,  Ks,  varies  from  0.  to  a maxi- 
mum value  cf  ± 1.75.  In  general,  it  is  sufficiently  accurate  to  round  off  the  *' 

value  of  Ks  to  the  fourth  decimal  place.  Therefore,  instead  of  using  a 

REAL*4  variable  to  store  KS,  an  integer  variable  of  2 bytes,  i.  e. , INTE- 
GERS, is  sufficient.  This  will  ameliorate  the  situation  of  storage  requirement. 

i 1 

c.  Storage  space  can  also  be  saved  if  one  makes  use  of  the  fact  for 
small  values  of  t,  the  table-size  is  smaller.  Hence,  a three-dimensional 
array  such  as  V{50, 50,40)  (see  Section  IV)  will  be  wasteful  in  storage  allo- 
cation. A better  alternative  will  be,  for  example  VI  (30,  30),  V2(32,  32), 

V3(36, 36). . . V40(50, 50).  This  would  amount  to  approximately  76,  000 
half-words,  or  38, 000  full  words  of  four  bytes  each. 
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*****  ••**»  PROGRAM  *TAS«S* 

PROGRAM  FOR  CALCULATING  THF  KERNEL  FUNCTION  KS(PX.  RY.  T> 

•••••  •••** 

COMMON/RNP ARM/ALPP . S7.0FLT * NT S * FLOPT (50*2) *«grp • EPo  c HT 
COMMQN/ST ATUS/ AT  ,BT . T I MF . ST  ATE ( SO . 3 ) * KERNEL (50*50) 

COMMQN/CNS  f AN/P I .PI 2*MACCY*  FACTOR 
COMMON/PR I PUN/ 1 R * I W *NOPRT 
COMMON  /SAVE/  R(4).  SI (4) * QF<4)*  NOROR(4) 

COMMON  /LIMIT/  AKU*  ALPHA*  3REAK*  AL2 
RE AL*8  PI*  024.  ANSI 20) 

REAL* 4 KERNEL .LENGTH 

DIMENSION  RX (50) « RY(30)*  TMS(30).  P(4*2) 

DIMENSION  SOLN(65)*  FRR(65>,  FSI (65) 

CALL  ERRSET  (208*999.-1*1) 

PI  = 3.1415926535398 
PI2  = .5*PI 
MACCY-6 
FACTOR  = 2. 

NOPRT  = -1 
ST-0.5 

ALPP  * .314159 
HT = 1 . 

ALPHA  = PI 2/ ALPP 
AL2  - ,25/(ALPHA**2> 

BREAK  =1.5 
XC  - 4. 

YC  = 2.5 
IACY  = 1 
ERR  = .05 

CALL  PREP 

VARIABLES  DEFINITION.. 

ST =BEAM/LENGTH  ratio 

ALPP  a PI/(2.*AIPHA) * WHERE  ALPHA  IS  NON-DIMENSIONAL  * I.E.  (LARRY 
DOCTOR'S  ALPHA) *A  * WHERE  A=  HALF-LENGTw  OF  CPAFT. 

RS*  alpha: *A 

HT  = DEPTH  OF  HALF-LENGTH  RATIO. 

ERR  * ERROR  ALLOWED  IN  THE  K- INTEGRAL 

TMS  = TIME  = NON-DIMENSIONAL  TIME  VALUE*  PRESENT = 0. 

= ACTUAL  TIME  INTERVAL  * S0RT(G/A) 

RX  = X-COORDINATE  OF  FIELD  POINT 
RY  a Y-COORDINATE  OF  FIELD  POINT 

IACY  IS  A PARAMETER  CONTROLLING  ACCUARCY  OF  THE  CALCULATIONS. 

I AC-Y=1  * FOR  ROUGH  CALC. 

IACY=2  FOR  GOOD  ACCURACY. 


READ  IN  ARRAY  OF  FIELD  POINTS 

READ  (5.1001)  NT,  (TMS(K),  K=1,NT) 
READ(S.lOOl)  NX.  (RXU),  I=1^NX) 
READ(5,1001)  NY.  (RY ( I ) ♦ 1=1, NY) 
READ  IN  ACCURACY  CONTROL  SIGNAL 
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READ  (5.  1033)  IACY 

1001  FORMAT  (I5/(F10.4)> 

1002  FORMAT  (F10. 4) 

1033  FORMAT  ( 15) 

IF  ( IACY  ,F0.  2 ) ERR  = .01 

* CALCULATE  APPROX.  UPPER  A<-  LIMIT  FOP  GIVEN  ERR. 

CALL  FINOKO  (AKU.  ERR.  ALPHA  ) 

00  100  K— 1 . NT 
TIME  * TMSfK) 

RD2  = HT*TIME**2  *i.l 
RD  = SORT (RD2) 


: PUNCH  OUT  GRID  SYSTEM  FOR  RECORD. 

WRITE  (6*1012)  TIME  • (RX(J).  J=l.  NX) 

1012  FORMAT  { lHO  • *■  NON-DIMENSIONAL  I ZED  TIME  VALUE  = •*  F10.5  • / 
1 IX.  *RX ( AT)  -».3X.  9F12.4  /(  12X.9F12.4)) 

WRITE  (7.  1013)  TIME 

1013  FORMAT  ( • KERNEL-TABLE  FOR  TIME  **♦  F10.5  ) 

WRITE  (7.1015)  NX.  ( RX(1).  1*1.  NX) 

1015  FORMAT  ( »RX  GRID.  NO.  OF  GRID  POINTS.  NX**  . 15  /{8F10.5)) 
WRITE  (7.1016)  NY.  ( PY ( J) • J=l.  NY) 


WRITE  (7.1016)  NY.  ( RY ( J) • J=l.  NY) 

1016  FORMAT  ( »RY  GRID.  NO.  OF  GRID  POINTS.  NY=» 
WRITE  (6.1024) 


15  /(8F10.5)) 


1024  FORMAT 


RY=BT*> 


DO  300  J*l.  NY 
BT  * RY(J) 

P ( 1 .2)  * BT.ST 
P(2.2)  * BT-ST 
R (3.2)  = PCI. 2) 
P(4.2)  = PC2.2) 
CALL  TIMING ( IJK) 


DO  200  1*1.  N) 

AT  * RX(I) 

PCI. 11  * AT*1 • 


P(2.1) 
P(3.X) 
P (4. 1 ) 


= P(l.l) 
= AT-1. 
= P (3. 1 ) 


SET-UP  CORNER  C00RI0INATF5. 

DO  120  M=1 .4 

R(M)  = SQRT (PCM.l )**2^P(M.2)**2) 

SI(M)  = 0. 

IF  ( ABS (P (M* 1 ) ) .GT.  1 .£-10  .OR.  ABS(P(M,2))  .GT.  l.E-10  ) 
1SI (M)  = ATAN2(P(M,2) . P(M.1>) 

120  CONTINUE 


FRR(I)  = SORT  (AT**2  ♦ BT**2) 

IF  { FRR ( I ) .LT.  0.5  ) FRR(I)  = 0.5 
FSI(I)  - 0. 

IF  ( ABS (BT)  .GT.  I.E-lO  .0».  ABS(AT) 
1FSKI)  = ATAN2(BT.  AT) 


.GT.  l.E-10  ) 


CHECK  IF  GRID-POINT  EXCEEDS  WAVE-FRONT. 

IF  ( BT  .LE.  YC  .AND.  AT  .GE.  (XC*RD))  GO  TO  3in 
IF  ( AT  .LE.  XC  .AND.  BT  .GE.  (YC*RD) ) GO  TO  3 IP 


IF  {((BT  .GF. 


•AND.  (AT  .GE.  XC)) 
98 
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C C CBT-YC) **2  ♦ (AT-*0>**2)  .GT.  «D2>)  60  TO  31 n 


CALCULATE  THE  INTEGRAL  9Y  CALLING  *1NTDK* 
VALUE  RETURNED  IN  ARRAY  *SOLN* 

CALL  INTDK? CFRR C T ) • SOLNCD*  UCY  ) 

GO  TO  200 

310  SOLNCI)  =0. 

200  KERNEL U*J)=  SOLNCI) 


CALL  TIMING  CIJD 
SECS  * #01*  C I JL-I JK) 

PRINT-OUT  LOOP. 


DO  400  M=1 « 10 

Ml  = CM— 1 ) *9*1 
M2  =M1*8 
IF  ( M2  #GT.  NX) 


M2=  NX 


IF  ( M #NE, 


1 1 GO  TO  420 

WRITE  C6*1003)«T*M,CS0LNCL)  • L-Ml,  M2*, 

1003  FORMAT  ClX*  F7.4*  16#  9F12.4  > 

WRITE  C7*1005)J*M*  (SOLNCL)#  L=M1.  M2) 

1005  FORMAT  ( 15#  13*  9F8.S) 

GO  TO  410 

* 

420  WRITE  (6*1006)M  (SOLNCL)*  L=M1*  M2) 

1006  FORMAT  CSX*  16*  9E12.4) 

WRITE  C7*1007)M*  C SOLNCL)  • L=M1*  M2) 

1007  FORMAT  C 18*  9F8.5) 

410  IF  C M2  #E0#  NX)  WRITE C6*  1004)  SECS 

1004  FORMAT  C 1H**  121X,  F8.2  ) 

IF  C M2  #GE.  NX)  GO  TO  300 

400  CONTINUE 

300  CONTINUE 
100  CONTINUE 
STOP 
END 


s ' 
! . 
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SUBROUTINE  !NTOK?(PR.  ANS*  KPT) 

REAL*8  VK • AKO.  ANSI*  ANS?*  AK 
REALMS  ANS3*  ANS4 

COMMON/RNPARM/ALPP « f T *DELT  ♦ NT  S . FLDPT  C 50  * ? ) , mFp  , poo  * h T 
REAL**  KERNEL 

COMMON/STATUS/AT. BT. TIME, STA1E (50*3) .KERNEL (50*50) 

REAL*8  PI 

COMMON/CNSTaN/  PI *PI2.«ACCY. FACTOR 
COMMON/PR I PUN/ I R • I W • NOP»T 
COMMON  /LIMIT/  AKU.  ALPHA*  BREAK*  AL2 
REAL *8  WW*?7*  W,  Z 

COMMON  /LAGUER/  Z?(68)«  WW(68).  Z(30>*  W(30)*  KHIGH.  Kl_OW 
DIMENSION  AUX (20) • KOUNT(IO) 

DIMENSION  ANSS (10) 

COMMON  /SAVE/RAO (4) . PST (4) *QF (4) «N3RDR(4) 

REAL*8  WNIMG 
EXTERNAL  WNING 

SUBROUTINE  FOR  CALCULATIMG( ALPHA) **-2  TIMES  THE  INTEGRAL  OF 
K=0.  TO  K=INFINITY  OF.. 

K*GAMMA*SIN(GAMMA*T)  *024  ( RIJ.  SITJ.  K)  * DK 
WHERE  GAMMA  = SORT (K*TANH(K*HT) ) • 

RR  a DISTANCE  OF  FIELD  POINT  FROM  CRAFT. S ORIGIN.*  TIME=T* 
AKU  SOHETIMES  REPLACES  THE  INFINITE  UPPER  LIMIT. 

ANS  - VALUE  OF  K-INTEGRAL*  R(I*J)  ANO  SKI.J)  IS  TRANSMITTED 
THROUGH  THE  COMMON  BLOCK  /SAVE/ 


REGION  OF  INTEGRATION  IS  DIVIOE0  INTO  TWO  PORTIONS* 
DIVISION  POINT  IS  GIVEN  BY  *8REAK* 

SET  *KPT*  EQUAL  TO  1 FOR  ROUGH  CALCULATIONS. 

SET  *KPT*  EQUAL  TO  ? FOR  REFINED  CALCULATIONS. 

N=1 

ANS  = 0. 

AKUP-O. 

PERIOD  = 2.*PI/ (PRESORT (HT)*TIME) 

GO  TO  (1.2).  KPT 

BLOCK  1 **** 

1 CONTINUE 

REGION  LESS  THAN  BREAK.  FREOUENCY=  R*5QRT (HT> *TIME 
REGION  BIGGER  THAN  BREAK.  FREQJENCY  = (R*TIM£/SGRT ( A*) ) 


00  150  L=1 *20 
AKLO=  AKUP 

AKUP  * AKLO  ♦ PERIOD 

IF  ( AKLO  *GE.  BREAK  I GO  TO  160 

IF  ( AKUP  .GT.  BREAK  ) AKUP=  BREAK*  .00001 

KGUNT ( I ) =6 

CALL  Q66  (AKLO*  AKUP,  WNING*  AN5S(1)) 
ANS  = ANS*  ANSS(l) 

IF  ( NOPRT  .N E.  1 ) GO  TO  150 
WRITE  (6*  1002)  L*  AKLO*  PERIOD*  ANSS Cl) 
1002  FORMAT  C 120*  2F10.4*  F12.6  ) 

150  CONTINUE 
GO  TO  161 


i 'it''.  :.rl  A 


c 
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160  AKU?=  BREAK 

161  3NS  - 0. 

DO  50  L=1 * 65 

AKLO  = AKUP  ' 

PERIOD^  2.*PI/CR9*TIMF/SQRT(AKLO> ) 

AKUPs  AKLO^PERIOO 

IF  ( AKLO  .GE.  AKU)  GO  TO  60 

IF  { AKUP  .GT.  A.KU  ) AKUP  = AKU*  .00001 

KOUNT { 1 } = 4 

CALL  OG4(AKLG*  AKUP*  WNTNG*  ANSS<2)) 

BNS  = BNS  ♦ ANSSt2) 

IF  ( NOPRT  .NE.  1)  GO  TO  50 
WRITE  (6*  1002)  L.  AKLO*  PERIOD*  ANSS(2) 
50  CONTINUE 

60  ANS  = AL2* (ANS*8NS) 


I 

c 

c 

RETURN 

I 

c 

BLOCK  2»**** 

2 CONTINUE 

00  350  L=1 *20 

1 

AKLO*  AKUP 

AKUP  s AKLO  ♦ PERIOD 

£ 

O 

IF  ( AKLO  .GE.  BREAK  ) 

IF  i AKUP  .GT.  BREAK  } 

KOUNT(l)  slO 

CALL  QG10  (AKLO*  AKUP 

; 

ANS  = ANS*  ANSS ( 1 ) 

IF  ( NOPRT  .NE.  1 ) GO 
WRITE  (6*  1002)  L*  AKLO 

350  CONTINUE 
60  TO  361 

360  AKUP  = BREAK 

361  BNS  =0. 

OC  450  L=1 • 65 
AKLO  = AKUP 

PERIOD*  2.*PI/CRR+TIHE/S0RT(AKL0)) 

AKUP=  AKLO*PERTOD 

IF  ( AKLO  .GE.  AKU)  GO  TO  460 

IF  ( AKUP  .GT*  AKU  ) AKUP  = AKU*  .00001 

KOUNT (11  = 7 

CALL  067 C AKLO*  AKUP*  WNTNG*  ANSS(2U 
BNS  s BNS  ♦ ANSS12) 


IF  ( NOPRT  .NE.  1)  GO  TO  450 
WRITE  (6*  1002)  L*  AKLO*  PERIOD*  ANSS<2) 
450  CONTINUE 
460  ANS  * AL2*(ANS*SNSr 


RETURN 

END 
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SUBROUTINE  FINOKO  (AKU.  ER.  ALPHA) 

OETERMINE  THE  APPROX.  UPPER  LIMIT  FOR  THE  INTEGRATION  WITH  RFSPECT  TO  K 
; VALUE  RETURNED  THRU  AKU. 

E = ER/ { SORT (ALPHA) *.5079491 ) 

AKU  = 5. 

OAK  =0.5 
DO  10  1=1.30 

IF  ((  SORT (AKU)  *EXP ( -AKU) ) .LE.  E)  GO  TO  30 
10  AKU  = AKODAK 
30  AKU  = AKU5*  .636620*  ALPHA 
: AKU  IS  K-SUB.  0. 

WRITE  (6.1001)  ER.  AKU 

1001  FORMAT  ( »0  ER=‘.  E10.3.  • AKU=* * F10.3  ) 

RETURN 

END 


U 


I 


I 


'/  4 
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FUNCTION  WNING ( AA  > 

EVALUATES  THE  INTEGRAND  OF  THE  K-INTEGRAL. 

NOTE..  THE  INPUT-  PARAMETERS  ARE..  AA*  TIME*  ALPHA.  HT * AND 

R<4)  AND  SKA)  IN  SAVE. 

REAL*4  KERNEL 

REAL*8  Q24  _ _ 

COMMON/RNP ARM/ALPP ♦ ST  * DELT ♦ NTS  * FLDPT (50*2)*  NFP  *ERR ♦ HT 
COMMON/ST ATUS/ AT  *8T  * T I ME  * ST  ATE (50*3)*  KERNEL (50*50) 

COMMON  /SAVE/RAD (4) * PSI (4) «QF(4) »NORDR(4) 

GAMMA  = SORT  <AA*TANH (AA*HT) ) 

WNING  = AA*GAMMA*  SIN<  GAMMA*TIME)  * Q24<AA*  ALPP) 

NOTE  THAT  Q24  TAKES  REAL*4  ARGUMENTS. 

RETURN 

END 


rrr.M  .. 


Computations  of  the  Kernel  Function 


Kernel 


Kernel-Table  for  TIME  * D. 50000 


Figure  B.4  KERNEL- TM3LE  FDR  TIME  «=  1 . SDOOO 


KERNEL ~7*BLE  FOR  TZME  * 5DDD0 


KEPNEL-’T^BLE  FDP  TIME  «=  3.00000 


KERNEL -T/iSLE  FDR  TZME  - 3. SDOOD 
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THE  THETA-INTEGRAL  ROUTINE  PACKAGE,  Q 24.  Documentation  is 
provided  here  for  a package  of  subroutines  that  calculates  the  theta-integral, 
Q(k,  R.'f').  The  mathematical  basis  behind  the  evaluation  of  this  integral 
has  already  been  given  in  the  second  section  of  this  appendix.  Consider 
equations  (136)  and  (145)  and  we  define: 


Q24(k,R  ,R  , 

* J 


a»)  = 2 C(k‘)[V  yj-  (kR.Jsin  (2n  % )], 

£=1.3.5... n M fl2*  13  13 


(153) 


where 


C (k‘)  = 
n 


T/2 


sin  2n0  d.  f 


sh(k*cos0)  sh(k*sin0) 


» k‘^0 


(154) 


R.. 

13 


ij 


{lKx-<-)lJ2  + [R ,.-(->V]2  }1/2 


X 


i » 1,2 
j = 1,2 


(155) 


(156) 


k« 


a‘ki 


(157) 


and  a'  * 7r/2a.  Equations  (153  - 157)  have  already  been  given  earlier. 

Here  then  are  repeated  for  the  sake  of  completeness.  Note  that  Q24  when 
multiplied  by  ky*  sinT  t,  as  can  be  seen  from  (136),  is  the  integrand  of 
the  k-integral.  The  physical  parameters  that  need  to  be  specified  before 
(153)  can  be  calculated,  are  a',  the  decay  factor,  and  /x,  the  beam  to  length 
ratio.  In  the  function  subprogram  Q24,  the  argument  list  consists  of  the 
value  k,  and  the  value  of  a’;  the  quantities  R..  and  f..  are  transferred  through 
a common  block.  This  and  other  informationlire  detailed  in  the  following 
paragraphs: 

a.  Limitations  of  Parameter  values. 

k*  = [.001,  125.3 
kR„  * [0.,  750.0] 

•*6 

Absolute  accuracy  is  approximately  10  . 

b.  Contents  of  Package  and  Brief  Description. 

Routine  Q24  (AK,  ALPP).  This  routine  takes  input  from  the  common 
block/ SAVE/ (see  user-instructions  section)  and  from  the  argument  list.  It 
serves  as  the  calling  program  for  routines  GCNK4  and  BESSEL  as  well  as  sums 
the  function  name,  Q24. 
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Routine  PREP,  This  short  routine  is  to  be  called  by  the  calling 
program  that  uses  Q24.  It  should  be  called  once  before  Q24  is  used.  Its 
purpose  is  simply  to  load  the  spline-coefficients  data  into  a common  block 
/ COEFG/  so  that  it  can  be  later  used  by  GCNK4. 

Routine  BESSEL. . An  efficient  routine  converted  from  the  MATH 
LIB  version.  It  generates  a series  of  Bessel  functions  for  a given  argument. 

Routine  SESSEL.  While  BESSEL  uses  double -precis  ion  arithmetic, 
routine  SESSEL  which  is  identical  to  BESSEL,  uses  single  precision  arithmetic. 

Routine  NBESSL.  A function  subprogram  called  by  Q24  for  estimating 
the  number  of  terms  required  in  the  Bessel  series.  Equation  (141)  is  used 
in  this  case. 

Routine  GCNK4.  This  routine  utilizes  the  input  data  in  common  block 
/ COEFG/  to  obtain  values  of  the  coefficients  Cn(k')  by  the  method  of  interpolation. 

One  Data  Deck  with  a header  card,  containing  values  of  spline  coefficients 
which  will  be  loaded  when  routine  PREP  is  called. 

c.  Instructions  for  Sample  Usage. 


C Declarative  and  other  statements  in  calling  program. 

REAL  *8  PI,  Q24 

COMMON/  CNSTAN/  PI , PI 2,  MACCY,  FACTOR 
COMMON/  SAVE/  R(4),  SI(4),  QF(4),  NORDR(4) 

DIMENSION  P(4,  2) 

PI  *=  3.14159265359 
P12  = ,5*PI 
MACCY  = 6 
FACTOR  = 2 

CALL  ERRSET  (208,999,-1,1) 

CALL  PREP 

C Generate  R..  and  'I'..  values  for  a given  (R  ,R  ) value  (ST=Beam/ 
1313  x y 

C length  ratio).  P(L, M)  below  are  the  cartesian  components  of  R..,  t..: 

l3 

P(l;  1)  = RX  + 1. 

P(l,  2)  = RY  + ST 
P(2, 1)  = P(l,l) 

P(2, 2)  = RY  - ST 
P(3, 1)  = RX  - 1. 

P(3, 2)  = P(l,  2) 

P(4, 1)  = P(3,  1) 

P(4,2)  = P(2, 2) 
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C 

© DO120M-1.4 

R(M)  « SQRT  (P(M1,  )**2  + P(M,2)**2) 

120SKM)  * ATAN2(P(M,2),P(M,1) ) 

C Define  AK  and  ALPP  before  calling,  and  then 

QVALUE  * Q24  (AK,ALPP) 

Note:  (1)  R-.  and  are  transmitted  through  the  COMMON  block/ SAVE/. 

(2)  QF(4)  and  NORDR(4)  are  values,  upon  exit,  corresponding  to 
(-)^Q(kR..,  k')  and  the  number  of  terms  used  in  the 

series,  respectively. 

d.  Computation  Time.  On  an  IMB  370  machine,  each  call  consumes 
approximately  . 01  sec  for  a value  of  kR  * 5 and  may  go  up  to  . 05  sec  for 
the  value  of  kR  * 150. 

It  is  worthwhile  to  point  out  that  another  example  of  the  application  of 
Q24  can  be  seen  in  the  program  TASKS  described  in  the  fourth  section  of 
O this  appendix.  In  that  case,  the  main  program  calls  PREP  and  sets  up  the 

quantities  R. . and  ..  The  function  Q24  is  then  called  by  the  routine  WNING. 
*3  y 

Listing  of  the  Q24  package  is  given  in  the  following  pages. 


■*> 
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DOUBLE  PRECISION  FUNCTION  024 ( T AK* ALPHA ) 

C ALPHA  HERE  IS  ALPHA-PRIME  IN  PROG.  REPORT. 

C 

C CALCULATE  THE  THETA-INTEGRAL  USING  THE  BESSEL  SERIES  EXPANSION 
C AND  SUMMING  THE  FOUR  SERIES*  EACH  CORRESPONDING  TO  A REF.  POINT 
C AT  THE  CORNER  OF  THE  CRAFT. 

C INPUT  VALUES  OF  THE  RADIAL  DISTANCE*  R*  AND  ANGLE*  SI*  COMES 
C THRU  THE  COMMON  BLOCK  /SAVE/ 

C 024  * SUM  I*  ( I—l *4)  0FC-1.)**I  * { SUM  N*  <N=ODD)  OF 
C -l.*J2N(AK*R(I))  * SIN(2N*SI(I))  * CM (AK* ALPHA)  ) 

C WHERE  CN(K)  IS  OBTAINED  BY  CALLING  THE  INTERPOLATION 
C ROUTINE  *GCNK4* 

C 

REAL*®  TSI*  TS*  TC*  TST • S4*  C4*  SUMO*  AKR*  UNITY (4) 

REAL*8  Y (2000) * PI*  ASIMPS.  AFILON 
rtTMPNCTftM  nilllklV(4Afl) 

COMMON  /SCRTCH/  Y«  NO*  NF«  ASIMPS*  AFILON*  ACUBNT * ATNH*  JJ 
REAL*B  JB(1031) * SIGNA 

EQUIVALENCE  (JB(1) *Yil) ) * (DUMMY (1 ) *Y 11032) ) 

COMMON  /CNSTAN/  PI*  PI2*  MACCY*  FACTOR 
COMMON  /RESULT/  CK(500) 

COMMON  /SAVE/  RR<4>*  PSI (4) * QF(4),  NORDR (4) 

DATA  UNITY ( 1 ) * UNITY<2)*  UNITY(3).  UNITY(4)  /4.D0*-4„D0*-4.D0* 

1 4*00/ 

FIND  MAX.  VALUE  OF  PR  IN  RR  ARRAY. 

AK  - TAK 

IF  ( AMALPHA  .LT.  .001  ) AK  * .001001/ALPHA 
IMKR*1 
AA*  RR(1) 

00  10  I»2*4 

IF  ( RRII)  .LE.  AA)  GO  TO  10 
!MKR*I 
AA*RR(I) 

10  CONTINUE 
AKR  = AA*AK 
PK  = AK*ALPHA 

C 

C CHECK  THE  MAX.  VALUE  OF  AKP. 

IF  ( AKR  .LT^  750.  ) GO  TO  11 
WRITE  (6*  1001)  AK.AA,  AKR 

1001  FORMAT  ( //*  5X*  * AK=»*  F10.4, «RMAX*» * F10.4*  » AKR**,  F10.4*  /* 
1 • AKR  TOO  LARGE  FOR  BESSEL  FUNCTION  ARGUMENT.*) 

STOP 

C 

C USE  MAX.  AKR  VALUE  AS  INPUT  TO  NBESSL  TO  DETERMINE  THE  ORDER 
C OF  BESSEL  FUNCTIONS  REQUIRED  (HENCE*  THE  NO.  OF  TERMS)  IN  THE 
C SERIES.  THIS  WILL  THEN  BE  USED  IN  THE  GENERATION  OF  THE 
C COEFFICIENT  ARRAY  CK. 

C 

C NBSL  IS  RETURNED  TO  BE  THE  ABSOLUTE  INDEX  OF  SERIES*  IT  IS  EVEN. 

C ABSOLUTE  INOEX  IS  THE  ORDER  OF  8ESSEL  FUNCTION  BEYOND  WHICH  TERMS 
C ARE  NEGLIGIBLE. 

11  NORDR (IMKR)  * NBESSL (AKR*  MACCY) 

NTRMS  * (NORDR (IMKR) *2) /4  ♦ 2 
CALL  GCNK4  (1*  NTRMS*  PK) 
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C 

C MOTE..  THE  COEFFICIENTS  *C**  ARE  GOOD  FOR  ALL  4 KR  VALUES. 

C 

024=0. DO 
C 

C L-INDEX  OF  THE  LOOP  BELOW  IS  AASSOCIATED  WITH  EACH  CORNER  OF  CRAFT. 
DO  20  L=l«4 

IF  ( ABS  (RP(L) ) .GT.  1.F-I0  ) GO  TO  21 
SUM0=0« 

GO  TO  35 
21  AKR=AK*RR(L> 

SI  = PSI(L) 

IF  ( L .NE.  IMKR)  NORDR(L)  = NBESSLUKR,  MACCY) 

C 

C GENERATE  VALUES  OF  BESSEL  FUNCTIONS  IN  SERIES. 

IF  ( AKR  -240.)  45.  46.  46 

45  NMAX  = 400 
SKR  = AKR 

CALL  SESSEL ( SKP . NMAX • JB • - 1 • • I ERROR « 1 • E-50 * DUMMY ) 

GO  TO  47 
C 

C SINGLE  PRECISION  BESSEL  FUNCTION  FOR  ARGUMENT  LESS  THAN  250.  . 

46  NMAX  = 1031 

CALL  BESSEL  (AKR.  NMAX.  JB.  -l.DQ.  IERROR.  l.E-50) 

C 

C SUM  BESSEL  SERIES 
C J(N)  IS  STOPEO  IN  JB(N*1) 

47  NN  = N0RDR(L)/2 
C 

TSI  = 2.D0*SI 
TS  * DSIN  (TSI) 

TC  * DCOS(  TSI) 

54  • 2.D0*TS*TC 

C4  = 1.D0  - 2.D0* (TS**2) 

SUMO  = JB (3)  * CK  ( 1 ) *TS 
J=1 
C 

C LOOP  INDEX  IN  SUMMING  OF  SERIES  IS  ODO  IN  VALUE.  STEP=2 
DO  30  N=3*  NN.  2 

C CALCULATE  SIN(2N*SI) 

TST  = TS 

TS  * TS*C4  ♦ TC*S4 
TC  * TC*C4  - TST*S4 
J = J*1 

30  SUMQ  = SUMQ  ♦ JB(2*NM)*TS*CK(J> 

C 

35  SUMQ  * SUMQ*UNITY(L) 

QF(L)  * SUMQ 
Q24  = Q24  ♦ SUMO 
20  CONTINUE 
C 

RETURN 

END 
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READS  IN  SPLINE  COEFFICIENTS  FROM  DATA  CARDS. 

COMMON  /COEFG/  NPK.  NPKMl * NCCAL.  PKSPClOO) «CSP(4.100.10) 

C 

READ  (5.1001)  NPK.  NCCAL 
NPKMl  * NPKsl 

READ(5. 1002)  (PKSP(K),  K=l*  NPK) 

READ  (5.1003)  (HCSP(J.K.N).  J*l*4),  K=l.  NPKMl).  N*l.  NCCAL) 

1001  FORMAT  (1615) 

1002  FORMAT  (6X.  6E11.4) 

1093  FORMAT  (6X.  4E16.7) 

RETURN 

END 
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FUNCTION  NBESSL  < AKR*  M) 

REAL*8  AKR 

C THIS  FUNCTION  DETERMINES  THE  NO.  OF  TERMS  REQUIRED  IN  THE  BESSL  SERIES  FO*j 
C GIVEN  ACCURACY  OF  10**<-M)  , AND  ARGUMENT 
C MUST  SOLVE  FOR  N SUCH  THAT  ( ALOGIO(N)  > ALOGIO(F)  ♦M/N 

C 

IF  I AKR  .GT.  l.E-5  ) GO  TO  10 

NBESSL =10 

RETURN 

10  F s 2.7182*AKR*.S 
ALF  = ALOGlOtn 
AM2  * M 

C 

NS»  1.36* AKR 

IF  ( NS.LT.  10)  MS*10 

DO  20  NsNSt  1000*  4 

IF^ALOGIO (FLOAT CN))*I1.*.5/N)  .GT.  (ALF*AM2/N))  GO  TO  30 
20  CONTINUE 

1001  FORMAT (t/*l*0IiATCH  - OUT  1*  REOUIRES  NBESSL  >1000*  BUT  SET=1000«) 

30  CONTINUE 

NSs  NBESSL/? 

IF  ( NS*2  .NE.  NBESSL)  NBESSL*NBESSL-1 

return 

END 
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SUBROUTINE  GCNK4  (Nl*  NTRMS*PPK) 

C THIS  SUBROUTINE  GENERATES  THE  ARRAY  CNK(PK)  AND  STORES  IT  IN  CK(N>. 

C HERE  N ACTUALLY  TAKES  ON  THE  VALUE  OF  N*1 *3*5*7  •••  ETC*  EVEN  THE 
C STORAGE  IS  IN  SEQUENCE. 

C VALUES  OF  CK( J) « J=N1 • NTRM,  ARE  CALCULATED  FROM  THE 
C EXISTING  SPLINE  COEFFICIENTS  STORED  IN  THE  ARRY  CSP. 

REAL*8  PI 

COMMON  /CNSTAN/PI*  PI2.  MACCY*  FACTOR 

COMMON  /COEFG/  NPK*  NPKMl*  NCCAL*  PKSP(IOO) *CSP(4.100.10) 

COMMON  /RESULT/  CK(SOO) 

C 

PK  * PPK 

IF  ( PK  .LT.  PKSPC1I)  GO  TO  100 
IF  ( PK  .LE.  PKSP(NPK) ) GO  TO  200 
100  CONTINUE 

C 100  WRITE  (6*10011  PK*  PKSP(NPK)*  ? KSP  (1) 

1001  FORMAT  (//*  • WARNING  REQUESTED  VALUES  OF  PK  IS  OUT  OF  THE  RANGE 
10F  PK  AVAILABLE*  ••  /*  • EXPRAPOLATIGN  IS  USED*  PK*»*  El 0.5* 

2 •PKSP(NPK)** *E10.5**PKSP < 1)  = • * E10.5  1 
C 

C SEARCH  FOR  LOCATION  OF  PK  IN  PKSP  ARRAY. 

200  DO  10  J=l*  NPKMl 
JJ  = J 

IF  C PK  .GT.  PKSP(J*1H  GO  TO  10 
GO  TO  20 
10  CONTINUE 
C 

• 20  NN  - NTRMS 

IF  I NN  .GT.  NCCAL 1 NN  = NCCAL 
HI  * PX-PKSP(JJ) 

H2  * Hl**2 
H3  = H2*H1 
PKSHP  = PK*SINH(PK1 
C 

00  30  N*N1*  NN 

C INTERPOLATE  CNK  FROM  TABLE  FOR  RANGE  OF  N SPECIFIED. 

CCC  = CSP(1*JJ*N)*H3  ♦ CSP(2.JJ*N)*H2  ♦ CSP(3*JJ*N)*Wl  * 

1 CSP(4.JJ*N> 

30  CKCN)  * EXP (CCC)/  PKSHP 


IF  ( NTRMS  .LT.  NCCAL  ) RETURN 
NN  * NCCAL* 1 
DO  40  N*NN.  NTRMS 

40  CK(N>  = PI2*DTANH(N*PI/PPK)  / PKSHP 
RETURN 
END 
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SUBROUTINE  BESSEL (X*N*APRAY*SIGNA*1ERR0R*  ALOW) 

PURPOSE 

TO  COMPUTE  ALL  ORDERS  OF  THE  BESSEL  FUNCTIONS  J(M.X)  OR 
E**(-X)*IfM*X)  FOR  A GIVEN  ARGUMENT  X. 


DESCRIPTION  OF  PARAMETERS 

X — THE  ARGUMENT 

N = THE  INTEGER**  MAXIMUM  DIMENSION  OF  ARRAY  IN  WHICH  THE 
RESULTS  ARE  TO  BE  PLACED  <0N  INPUT) 

N st  THE  INTEGER**  INDEX  OF  THE  LAST  ELEMENT  OF  ARRAY  WHICH 
CONTAINS  A NON-ZERO  RESULT.  (ON  OUTPUT) 

ARRAY  s THE  REAL*8  ONE-DIMENSIONAL  ARRAY  OF  LENGTH  N IN 
WHICH  THE  RESULTS  APE  TO  BE  PLACED.  AFTER 
EXECUTION  OP  THIS  SUBPROGRAM.  THE  MTH  ELEMENT  OF 
ARRAY  CONTAINS  JfM-i.X)  OR  E**(-X)*IfM-l,X) . 

SIGNA  as  THE  REAL*8  INDICATION  SPECIFYING  WHICH  KTN0  OF 
BESSEl.  FUNCTIONS  APE  TO  BE  COMPUTED.  SET  SIGNA 
= -I  FOR  -JCM.X)  OR  SIGNA  = 1 FOR  E**t-X)*I fM.XI . 
IERROR=  AN  ERROR  INDICATOR  RETURNED 
1 ERROR  * 0 MEANS  NO  ERROR 

I ERROR  * 1 MEANS  ARGUMENT  TOO  SMALL*  ABSfX)  .LE.  10**-38 
I ERROR  ” 2 MEANS  ARGUMENT  TOO  LARGE*  ABSfX)  .GT.  ?S0 
FOR  J(M*X)  OR  ABSfX)  .GT.  350  FOR  E**f-X)*HM.X) 

ALOW  * REAL**  VARIABLE  FOR  ERROR  CONTROL.  FOR  BEST 
ACCURACY.  USE  l.E-78. 

DESCRIPTION  OF  PROGRAM 

WHEN  ABSfX)  .LT.  0.1*  THE  BESSEL  FUNCTIONS  ARE  COMPUTED 
USING  THEIR  POWER  SERIES  EXPANSIONS. (SEE  THE  HANDBOOK  OF 
MATHEMATICAL  FUNCTIONS.) 

WHEN  ABSfX)  .GE.  0.1.  THE  BESSEL  FUNCTIONS  ARE  COMPUTED 
USING  A RECURSION  RELATION  METHOD.  CSEE  “GENERATION  OF 
BESSEL  FUNCTIONS  ON  HIGH  SPEED  COMPUTERS") 

REAL*8  ARRAY (1).  ZUM*  SIGNA*  X 

OATA  F2/1. 3591*1/. 8L0W/1.E-38/ 

DATA  IFLG/0/.  NMAX/1030/. BREAK/. 1/ 

IF(IFLG.NE.O)  GO  TO  98 
DLOW*ALOG< l ,E14*AL0W) 

IFL6*l 
IERRORsO 
NM1  *N-1 
SIGNTsl. 

IFfX.LT.O)  SIGNT=-SIGNT 
X*DABSfX) 

XX  * X 

DO  99  I*I*N 

ARRAY ( I ) *0.09 

IF (XX.GE. BREAK)  GO  TO  200 


SERIES  METHOD  - USED  FOR  X SMALLER  THAN  .1 

IF  ( XX  .GT.BLOW)  GO  TO  101 

ARRAYU)  = 1.00 

N=1 
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IERR0P=1 

RETURN 

101  CLOW=10.*ALOW/X 
EXPX=1. 

IF (SIGNA  .GT.0.D0)EXPX=EXP(XX) 

1=1 

FMULT=X*X*SlGNA/4. 

E=1  • 

GO  TO  103 

102  E=E* (X/2. ) /FLOAT  ll-l) 

103  SUM=1. 

ANPLt'SsFLOAT  ( I ) 

AN=1. 

TERMsSUM 

104  SUME=SUM 

TERM=TERM*FMULT/ ( ANPLUS*AN) 

SUMsSUN^TERM 

ANPLUS=ANPLUS*1. 

AN=AN«1. 

IFCABStSUME-SUM)  .GT.  l.E-70)  GO  TO  104 
ARRAY  C ! ) =E*SUH/EXPX 
1*1*1 

IF  (I  .LE.  N .AND.E  .GT.  CLOW)  GO  TO  102 
N=I-1 
GO  TO  205 
150  IERR0R=2 
N*0 

RETURN 

C 

C RECURSION  RELATION  MEHTOD  - USED  FOR*  X LARGER  THAN  .1 
C 

200  CLXE  = ALOG(XX*F2l 

IF  CSIGNA  .LT.  0.D0)  GO  TO  207 

IF  (XX.GT.  750.)  GO  TO  150 

CLOW=DLOW*X 

I=IFIXCXXi 

GO  TO  201 

207  IF  (XX.GT • 750.)  GO  TO  150 
ClOW=DLOM 
I = IFIXC1.38*XX) 

201  1*1*3 
AN=FLOAT(I) 

IF  ( (-ALOGCAN)*t.5»AN)+AN*CLXE)  .GT.  CLOW  .ANO.  U.LT.NMl)  >GOTO  201 

IF (I  .GT.NMAX)  GO  TO  150 

N=I 

CC  SS=-AL0GfAN)*(.5*AN)  ♦ AM*CLXE 

C WRITE (6.  1001)  SStCLOW.  CLXE 

C1001  FORMAT  ( l HO  * * CHECK-  »,3E15.6) 
array  < i ) =1 ,n-64 

00  202 -1=2. N 
J=N- 1*1 

202  ARRAY (J)=ARP AY ( J*2)*SIGNA*2.DO*OFLOAT  C J) *ARRAY  { JM ) /X 
2UH=ARRAY (1) 

J=1 

IP (SIGNA  .UT.0.00) J=2 
JLOW=J» 1 

DO  203  I=JLOW*N* J 
203  ZUM=  ZUM*2.nO*ARRAY(I) 

DO  204  T=l*H 

204  ARRAY ( I } =APRAY ( I ) /ZUM 
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?05  IF(SIG*T  ,GT.O»  RETURN 
00  206  I=2.M *2 
?06  ARPAY(I)=-AR»AV(I) 
RETURN 
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SUBROUTINE  SESSEL  C X * N * ARRAY • S I SNA • I ERROR*  ALOW*  SPRAY) 

C 

C ROUTINE  SIMILAR  TO  *8ESSEL**  EXCEPT.. 

C SINGLE  PRECISION  ARITHMETIC  IS  USED. 

C *BRRAY*  IS  USED  IN  COMPUTATIONS.  * ARRAY*  RETURNS  THE  RESULTS. 

C NOTE.  SIGMA  IS  REAL *4*  X IS  ALSO  REAL***  WHILE  * ARRAY*  IS  REAL *8 

C 

REAL*8  ARRAY (1) 

DIMENSION  RRRAYC1) 

DATA  F2/1.359141/.BL0W/1.E-38/ 

DATA  IFLG/©/*  NHAX/  400/*8REAK/.l/ 

IFIIFLG.NE.O)  GO  TO  98 
DL0W*ALO5{ 1 .fel4*ALOW) 

IFLG*1 

98  I ERROR* 0 - 

NMl  *N-1 
SIGNT*1. 

IF(X.LT.O)  SIGNT=-SIGNT 
X*  ABS(X) 

XX  = x 
DO  99  1=1,  N 
99  BRRAYIDxO. 

IF (XX.GE.BREAK)  GO  TO  200 
C 

C SERIES  METHOD  - U5E0  FOR  X SMALLER  THAN  .1 
C 

IF  ( XX  .6T.BL0W)  GO  TO  101 

ARRAY Cl 1*1. DO 

N=1 

IERR0R=1 

RETURN 

101  CL0N=10.*AL0W/X 
EXPX*1. 

IF (SIGN*  .GT.O.  )EXPXsFXPCXX) 

1=1 

FMULT=X*X*SI GNA/4 . 

E*l. 

GO  TO  103 

102  E=E*CX/2.) /FLOAT (1-1) 

103  SUM*1. 

ANPLUS*FLOAT  { I ) 

AN=1  • 

TERM=SUM 

104  SUME=SUM 

TERM=TERM*FMULT/ ( ANPLUS*AN) 

SUM=SUH*TERM 

ANPLUS=ANPLUS*1. 

AN=AN*1. 

IFCABS(SUME-SUN)  .GT.  l.E-70)  GO  TO  104 
ARRAY  C I ) =E*SUM/EXPX 
1*1*1 

IF (I  .LE.  N .AND.E  .GT.  CLOW)  GO  TO  102 
N*I-1 
GO  TO  205 
150  IERR0R*2 

WRITE  (6*1001)  X*N 
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1001  FORMAT  { • STOPPED  IN  SESSEL.  X*,  F10.4,  » NREOO  =*.12. 

1 * .GT.  NMAX*  ) 

N=0 

STOP 

RECURSION  RELATION  HEHTOO  - USED  FOR  X LARGER  THAN  .1 

200  CLXE  = ALOG(XX*F2) 

IF  CSIGNA  .LT.  0.  ) GO  TO  207. 

IF  (XX. GT.  350.)  GO  TO  150 

CLOW=DLO¥*X 

I=IFIX(XX) 

GO  TO  201 

207  IF  (XX.GT.  ?50.)  GO  TO  150 
CLOV-DLOtf 
I = IFIXCl .?8*XX) 

201  1=1*3 
AN=FLOAT { I ) 

IF  ((-AL0G(AM)»(.5*AN)*AN«CLXE>  .GT.  CLO*  .AND. (I .LT.NMl )) GOTO  201 

IF  (I  .GT.fiU&X)  SO  TO  150 

M=I 

BRRAYCI)  = 1.E-64 
BRRAY(I*1)  = 0. 

DO  2C2  1=2. N 
J“N— 1*1 

20 2 8RRAY (J) =SRRAY I J*2) *SI6NA*  2.*  FLOAT CJ) *8RPAY (J* 1 )/X 
ZUM  = 3RRAY ( 1 ) 

J=1 

IFCSIGNA  .LT.O.  ) J=2 
JLOM=J*l 

DO  203  I=JLOtf«N*J 
203  ZU*=  ZUK*2.  *BRRAY(I) 

DO  204  1=1 .X 

204  ARRAY { I )=BROAY< I) /ZUM 
?05  1F(S1GNT  .GT.O)  RETURN 
DC  206  1=2. N, 2 
206  ARRAYCI)s-BRRAY(I) 
return 
END 
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